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Abstract 


We  categorize  some  of  the  finite  difference  methods 
that  caji  be  used  to  treat  the  initial  value  problem  for 
the  boundary  layer  differential  equation  ixy'  =  f(y,x). 
These  methods  take  the  form 


k  k 

—  '"i  ^n+i  =  ^^"^  ^ 
1=0   ^  ^^^  i=0 


C  '^i  ^n+i   =  ^^"''C  ^i  ^^  Vi'^n+i^' 


where  ix  =  h^,   0  <  7  <  1,   and  h  is  the  meshsize.   We 
define  a  new  kind  of  stability  called  M--stability  and  prove 
that  under  certain  conditions  tx-stability  implies  convergence 
of  the  difference  method.   We  investigate  ii-stability  and 
the  optimal  methods  which  it  allows  (i.e.  methods  of  maximum 
accuracy) .   These  results  are  generalized  to  the  vector 
case. 

We  also  consider  the  related  finite  difference  methods 
for  the  initial  value  problem 

l^y'   =  f{y,z,x) 
z'   =  g(y*z,x). 

To  illustrate  the  theoretical  results,  a  niimber  of  problems 
were-  solved  by  these  methods  with  the  use  of  New  York  Uni- 
versity's CDC66OO  and  IBM709^  computers. 
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I.      NUMERICAL  ANALYSIS  AND  BOUNDARY  LAYER  PROBLEMS 


Numerical  methods  for  solving  the  Initial  value 
problem  on  the  Interval   [0,1]   for  the  ordinary  differ- 
ential equation  (ODE)  problem 

(1.1)  y'   =  f(y,x)   and  y(0)   =  1 

are  well  understood.   Among  the  most  useful  methods  are 

the  finite  difference  schemes.   For  these,  we  divide  the 

interval   [0,1]   into  mesh  points  at  a  distance   h   apart 

X  =  0,    x-.    -  h,    Xr.   =   2h,  ....  X      =   nh.   An  important  class 
ol^  n 

of  finite  difference  methods  can  be  written  as 

k  k 

a.  Y  ^.   =  h  / 

i^  ^  ^+^       i=0 


(1.2)        }        a.  Y  ^.   =  h  /   p.  f(Y  _^.,x  _,.  ) 


where   k   is  a  fixed  integer  and  a   and   P    (v  =  0, 
l,...,k)   denote  real  constants  which  do  not  depend  upon 
n.   The  first  proof  that  the  solution  of  equation  (1.2) 
converges  to  the  solution  of  equation  (l.l)  if  and  only  if 
certain  algebraic  conditions  are  satisfied  was  given  by 
Dahlquist  [ 1] .   Additional  literature  may  be  found  in 
Henrici  [2]  and  Hull  and  Luxemburg  [3].   We  will  describe 
this  theory  briefly  in  order  to  prepare  the  reader  for  the 


more  complicated  nximerical  analysis  of  boundary  layer 
equations. 

We  associate  with  the  difference  equation  (1.2)  two 
polynomials 


p(0   =  a^C^  +  CL^_^^^-^  +   '••    +  a^ 


d(0    =   pj^c^  +  \.^^^~^  +  •••  +  Pq- 

We  assTime  for  convenience  that  p(C)   and  (5(C)   have 
no  common  factors.   In  order  that  the  difference  equation 
be  a  good  approximation  to  the  ODE  y'  =  f(y,x)   we  must 
prescribe  certain  consistency  conditions  on  the  coefficients 
a^     and     P^. 

DEFINITION.   A  difference  equation  of  the  form  (1.2) 
is  said  to  be  consistent  if 

a)  p(l)   =  0 
and 

b)  p'(l)   =  S(l). 

We  must  also  ensure  that  the  difference  equation  be 
insensitive  to  small  changes  in  the  local  errors.   This 
will  be  guaranteed,  as  shown  in  [ 1] ,  by  the  stability 
condition. 


DEFINITION.   A  difference  equation  of  the  form  (1.2) 
is  said  to  be  stable  if  the  roots  of  p(0  =  0  lie  in  or 
on  the  unit  circle  in  the  complex   C -plane,  and  are  simple 
if  they  lie  on  the  circle. 

Since  (1.2)  is  a  k    order  difference  scheme,  we 
must  supply  k  initial  values  Y  ,Y-,,...,Y,  ,   in  order 
to  solve  the  equation.   We  may  let  Y  =  y(0)  =  1   and 
calculate  by  Taylor  series,  Runge-Kutta,  or  some  other 
method  the  values  of  Y-, ,  Yp,  . . ., Y,  -,  .   We  merely  prescribe 
that  the  method  chosen  to  calculate  the  initial  data  allows 
us  to  write 


(1.5)   |Y^-y(x^)|   -   (^(h)   as  h-^0  for  v  =  0, 1,  . .  .,k-l. 


The  more  general  difference  equation 


k 


(1.4)      >   a.  Y  ^.   =  h  >   p.  f(Y  ^.,x  ^.  )  +  R 
^         4 — TT-  1  n+1       4 — 7T     1    n+i'  n+i     n 
1=0  1=0 


may  also  be  considered. 

Here,   R^   is  the  round-off  error  and  arises  from  the 
inaccuracy  of  the  arithmetic  process  used  to  solve  the  differ- 
ence equation.   To  guarantee  convergence  when  using  equation 
(l.'^)  we  add  the  condition  that 


^n  "  Oi^^)      as  h-^0, 


With  the  difference  equations  (1.2)  or  (1.^)  we  asso- 
ciate the  difference  operator 


(1.5)   L[y{x);h]   =  /__  [oi.    y(x+ih)  -hP  y  (x+ih)] 

i-0 


We  will  assinne  the  function  y(x)  has  continuous  deriva- 
tives of  sufficiently  high  order  and  expand  y(x+ih)  and 
y'(x+ih)   in  Taylor  series: 


y(x+ih)   =  y(x)  +  ihy'(x)  +:^i^h^y"(x) 


hy'(x+ih)   =  hy'(x)  +ih^y"(x)  + 


Substituting  these  expressions  into  equation  (I.5)  gives 


L[y(x);h]   =  C^y(x)  +C^hy'(x)  +C2h^y"(x) 


+  C  h^y^'^^x)  + 
q  ^ 


where 


C 


1 


-  (a   +  2'^ap  +  ...  +  k^^) 


-  O^  +  2^-^P2  +  ...  -^  k^"^?i,) 


(q-l)I 


We  note  that  by  the  consistency  condition  C   =  p(l)  =0 

and   C,  =  p '  ( l)  -  (3(  l)  =  0.   A  difference  operator  of  the 

form  (1.5)  is  said  to  be  of  order   p   if  C   =  C   =  •••  = 

C   =0,   but   C  ,-,  i^  0.   In  this  case  we  may  write 
p    '  p+1 


L[y(x);h]   =   (5-(hP"^^)   as  h  -^  0 

if  the  derivatives  are  also  bounded.  The  associated  differ- 
ence equation  satisfied  by  the  solution  of  (l.l)  is 

k  k 

)    a.  y  , .   =  h  /    p.  f  ^.  +  T 
^— pr  1  "^n+i       ^. — ^  1  n+i    n 
1=0  1=0 

where   T   is  the  truncation  error   (T  =    0  i^        )   3'S 

h-0),  y^+i  =y(x^+i),  and  Vi  =  f(y(x^^i),x^^.). 

We  are  now  in  a  position  to  define  convergence. 

DEFINITION.   A  difference  equation  of  the  form  (1.2) 
or  (1.^)  is  said  to  be  convergent  if 


lim  Y   =  y(x  )    for  all   x  e  [0,1]. 
h  -►  0  ^        ^ 
x^=nh 

The  main  theorem,  first  proved  by  Dahlquist  [ 1] ,  is 

THEOREM  1.1.   A  necessary  and  sufficient  condition 
for  convergence  of  the  consistent  finite  difference  equation 
given  by  (1.4)  is  stability. 

Proof:   See  [ 1] ,  [2] ,  or  [3] • 


Remark  1.  We  have  ass\amed  that  at  least  y"(x)  is 
continuous  and  bounded.  While  less  stringent  conditions 
may  be  used  for  the  equation  (l.l),  this  will  not  be  the 
case  when  we  study  the  boundary  layer  equation. 

Remark  2.  Although  we  only  considered  the  scalar 
equation  (l.l),  the  same  theory  is  applicable  to  the  vector 
equation  y'  =  f(y,x)   where  y'   and  f  are   s-dimen- 
sional  vectors. 

Boundary  layer  phenomena  may  be  associated  with  the 
equation 

(1.6)  M-y'   =  f(y,x)    and  ,  y(0)   =  1  for  x  e  [0,1], 

where  m-  is  a  small  parameter.  A  number  of  questions 
can  be  posed  for  such  a  system.   For  example, 

a)  What  is  the  nature  of  such  a  solution  for 
small  M-? 

b)  As  \i  -*■  0     does  the  solution  y(x,M')   to 
equation  (1.6)  approach  a  solution  to  the 
reduced  equation 

(1.7)  f(y,x)   =  0, 

and  if  so,  which  one? 

c)  Can  a  \miform  asymptotic  expansion  be  foimd 
for  the  solution  of  (1.6)? 
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These  questions  have  all  been  answered  in  recent  years  not 
only  for  the  equation  (1.6)  but  also  for  the  vector  system 


H  g  =  f(y,z,x)    y(0)  =  y° 


(1.7)         H  =  g(y,z,x)    z(0)  =  z° 


0  <_  X  <   1 

where  y  is  an  M-dimensional  vector  and  z  is  an  m- 
dimensional  vector.  We  will  briefly  describe  the  treat- 
ment found  in  Vasil'eva  [^] . 

We  first  introduce  some  definitions.   Let  y  =  (t)(z,x) 
be  one  of  the  solutions  of  the  system  f(y,z,x)  =0.  By 
the  degenerate  system  of  equations  corresponding  to  the 
root  y  =  (t)(z,x)  we  will  mean 

y  =  <t)(  z,x)      0  <_  X  j<  1 
(1.8) 

If  =  g(y,z,x)    z(0)  =  z°  . 

Let  the  solution  of  this  system  be  denoted  by  'z(x)   and 
y(x)  =  4)(z(x),x) . 

DEFINITION.   The  root  y  =  <|)(z,x)   is  isolated  on  a 
bounded  closed  set  D  =  ((z,x)}  if  there'  exists  an  e  >  0 
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such  that  the  system  f(y,z,x)  =0  has  no  solution  other 

that  (j)(z,x)   for   |y  -  (t)(z,x)|  <  e. 

o 
Example  1.   Suppose   f(y, z,x)  =  -y  +1.   The  roots 

are   y  =  ±1   and  are  obviously  isolated. 

2     2 
Example  2.   Suppose   f(yjZ,x)  =  -y  +  x  .   The  roots 

are   y  =  ±x.   These  roots  are  not  isolated  at   x  =  0   if 

we  are  considering  the  interval  [-1,1]   but  are  isolated 

on  the  interval  [-^,1]  . 

The  equations 


(1.9)  1^  -  f(y,z*,x*) 


( z*   and  X*   are  regarded  as  parameters)  will  be  called 
the  adjoined  system  of  equations. 

DEFINITION.  The  isolated  root  y  =  (t)(z,x)  will  be 
called  positively  stable  in  D  if  the  real  parts  of  the 
roots  of  the  characteristic  equation 


j^^^(||(<t)(z,x),z,x)  _  ^^)   ^   Q 


(  -v-  is  a  matrix)  are  negative  in  D. 

DEFINITION.  The  domain  of  influence  of  an  isolated 
positively  stable  root  y  =  ^{z,x)  is  the  set  of  points 
(y*,z*,x*)   such  that  the  solution  to  the  adjoined  system 


(1.9)  satisfying  the  initial  conditions   y   ^q  =  y*   tends 
to  the  value  ^{z*  ,x*  ) ,      as  t  -»■  oo. 

We  can  now  state  the  main  theorem  about  boundary  layer 
equations: 

THEOREM  1.2.   If  some  root  y  =  'Ii(",x)   of  the  system 

f(y, z,x)  =0   Is  an  isolated  positively  stable  root  in  some 

bounded  closed  domain  D,   if  the  initial  point   (y°,z°,0) 

belongs  to  the  domain  of  influence  of  this  root,  and  if  the 

solution   z  =  'z(x)   of  the  degenerate  system  (1.8)  belongs 

to  D   for   0  <  X  <  X,   then  the  solution   (^S^'l^!)   of  the 
—   —  ^z(x,M.) 

original  system  (l.?)  tends  to  the  solution   (^ U(x),x))   ^^ 
the  degenerate  system,  as  m-  "*"  0,   the  passage  to  the  limit 


lim  y(x,M.)   =  y(x)   =  ^{z{x),x) 
M-  -*  0 


holding  for  0  <  x  <_  1  <  X,   and  the  passage  to  the  limit 


lim  z(x,M.)   =  z(x) 
M-  -►  0 


for   0  ^  x  ^  1  <  X. 

Proof:   See  Vasil'eva  [^] . 

We  illustrate  the  theorem  by  considering  a  few  examples: 

Example  '^ , 

M-y'  =  -y   y(o)  =  i. 
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lated.   The  solution  is  y  =  e  '     . 


The  only  root  is  y  =  0  and  therefore  it  is  iso- 

The  graph  of  this 
solution  is  given  in  Figure  1.   The  dotted  curves  show 

the  solution  for  values  of 
M.-^0.   The  small  interval 
in  which  the  slope  of  the 
curve  is  changing  most  rapidly 
is  called  the  "boundary  layer. 
A  conservatives  estimate  of 
this  interval  is   [  0, -A|-Lin  |x] 
where   A   is  a  positive  con- 
stant that  is  independent  of 
M- .   For  the  scalar  equation 
l-'-y'  ~   f(yjx)   it  can  te  shown 


\ 


^ 


\^ 


1 

Fig.  1 


that  A  =  2/ L  if   -L  <_  f  (y,x)  <_  -L  <  0.   See  Vasil'eva 
for  details. 


Example  \ , 


M-y 


=  -y(y-l)(2x  + 1)10 


with  y{0)  =  2  over  the  interval   [0,1]. 

The  roots  of   f(y,x)  =  0   are  y  =  ^{x)    =0   and 
y  =  (t)(x)  =  1.   They  are  clearly  isolated.   The  initial  point 
(2,0)   belongs  to  the  domain  of  influence  of  the  root 
(t)(x)  =  1  because  the  adjoined  equation 
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^    =  -y(y  -  1)10 


with  initial  condition  y   _^  =  2   has  as  solution  y(T)  - 
2/(2-e     ).   We  see  at  once  that   lim  y(T)  -   1.   Purther- 

T  — »-oo 

more,    the   root      y  =  it)(x)    =  1      is   positively   stable   in      [0,1] 
"because 


^(l,x)       =      (-2y+l)(2x+l)l0      ^^=      -10(2x+l)       <      0, 


Therefore,  we  would  expect  that  the  solution  y(x,M.)— ^1 
as  M-~^0.   Indeed  this  is  the  case  because 


y^""'^^   =  2  _  g-10x(x+l)/.L 


and 


^^%    ,  _  ^-10x(x+l)/,  -     ^         if  X  ^  0, 


The  paper  by  Vasil'eva  [t-]  goes  en  to  explain  how  to 
find  an  asymptotic  expansion  of  the  solution  to  the  system 
(1.7)  in  terms  of  the  small  parameter  \x .      Her  in  addition 
to  the  conditions  of  Theorem  1.2  we  assume  tnat  the  righthand 
sides  of  (1.7)  have  continuous  partial  derivatives  of  order 
up  to  n+2.   With  this  condition  Vasil'eva  finds  an  asymptotic 
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expansion  for  0  Jl  ^  Jl  1  which  contains  n  terms.   Inside 
the  bo\indary  layer   [  0, -Aia-^nM.]  each  terra  contains  three 
functions  found  by  solving  three  separate  differential 
equations  or  transcendental  equations.   Outside  the  boim- 
dary  layer   (-AM.^nM-,1]  the  terms  are  much  simpler  and  can 
be  astermined  from  the  variational  equations.   This  proce- 
dure for  finding  the  asymptotic  expansion  is  a  very  tedious 
one  and  can  only  be  explicitly  calculated  for  the  simplest 
problems . 

It  is  the  aim  of  this  paper  to  tie  together  the  known 
numerical  analysis  theory  with  the  boundary  layer  theory 
in  such  a  way  that  this  problem  can  be  solved  with  computers 
even  as  \i-*-0  .      We  do  this  by  allowing  m-   to  vary  with  the 
meshsize  h.   That  is,  we  let  (j.  =  h^  where  0  <  7  <  1. 
(Of  course  for  each  fixed  value  of  \i,      the  Dahlquist 
theory  implies  convergence  as   h— »-0.   But  the  convergence 
is  not  uniform  in  \x,      since  the  derivatives  of  the  solution 
y{x,M.)   may  become  unbounded  as  m.-^0.) 

If  we  attempt  to  apply  the  proof  of  Theorem  1.1  to 
this  new  problem,  we  run  into  serious  difficulties  because 
the  proof  uses  the  fact  that 

Mx 

(1.10)  lira  (1  +  Mh)^   =  e   ^ 

h  -►  0 

X  =nh 
n 

where  M  is  a  positive  constant.   In  our  case,  however. 
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the  corresponding  quantity  that  would  ha.ve  had  to  be 
estimated  is 

1  _  -i    X  /h 

(1.11)    lim  (1  +  Mh  ^)^  =  lim  (1  +  Mh  ^^)  ^    ^     ^ 
h  ->■  n  h-^-O 

n 

as  can  be  shown  by  a  simple  application  of  L' Hospital's 
rule '      Let 

T    X  /h 
A   -  lim  (1  +  Mh-^"^)  "   . 
h^O 


Then 


^n  A  =     lim 


X   in(l  -r  Mh-^"^) 


h-.0        ^ 


X  (1  -  7)h"^  M 

=  lim  =j where   x  e  (0,1] 

h  — 0   1  +  Mh  ^^ 


=     +00 


or 


A   =  ot. 


instead  we  will  take  advantage  of  the  positive  stabi- 
lity property  and  make  use  of  the  fact  that 


■-,  -,    X  /■  n 

(1.12)   lim  (1  -  Mh"^"^)^  =  lim  (l  -  Mh"""/)  ^^    =  0. 

h  -^  0  h  -<■  0 

X  =nh  X  =nh 

n  n 


ih 


We  pay  a  price  for  the  privilege  of  using  (1.12)  instead 
of  (I.IO),  namely,  we  must  restrict  ourselves  to  a  smaller 
class  of  polynomials   6(C),   to  guarantee  convergence. 
This  class  will  contain  optimal  method  schemes  of  order  at 
most  k  -I  1   in  stead  of  k  -^  P.   as  is  the  case  with  the 
ODE  troate^l  "iv    .  ...^aist's  work. 

We  will  also  use  a  slightly  different  definition  of 
convergence  for  boundary  layer  difference  equations. 

DEFINITION.   A  difference  equation  associated  with  a 
"boundary  layer  differential  equation  will  "be  called  conver- 
gent if  for  all  h  sufficiently  small  (i.e.   h  <_  h^)   there 


exist  constants   T  and  7  <  1  such  that 


:Y   -  y  1   <   Th-'-"^ 


for  all  nh  =  x^  e  [0,1] . 

We  can  look  at  this  idea  of  relating  M-   to  h   as 
follows:   Given  a  differential  equation  M-y'  =  f(y,x)   we 
ask  how  can  we  choose   h  so  that  the  associated  difference 

equation  will  give  an  accurat  approximation.   If  p.   is 

1/7 
sufficiently  small,  we  choose   h  by  the  formula  h   =  [i  '  ' 

where  0  <  7  <  1.   Furthermore,  we  ask  what  is  the  nature 

of  the  solution  to  this  difference  equation  as  I-l— ►O.   (It 

is  obvious  that  if  M-~*"0  we  must  choose  smaller  and  smaller 

meshwidths  in  order  to  obtain  a  good  approximation  to  the 

solution  of  the  differential  equation.  ) 
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II.  LINEAR  EQUATIONS  AND  |a-STABILITY 


We  first  study  the  numerical  methods  associated 
with  the  linear  equation 

(2.1)  [ly'      =     -L^y  +  m(x) 


on 


the  interval   [0,1]   with  initial  conditions   y(0)  =  1. 


L   is  a  constant  which  is  greater  than  zero  and  m'(x)   is 
a  continuous  function  of  x.   As  we  shall  see  the  analysis 
of  convergence  will  be  independent  of  m(x)   because  we 
are  interested  in  only  the  difference   i Y^  -  y(x  ) | .   We 
therefore  consider  the  homogeneous  equation 


ixy'   =  -L  y  for   0  <_  x  <_  1 
(2.2) 

with  initial  conditions   y(0)  =  1. 

It  will  turn  out  that  the  results  obtained  for  the 
linear  equation  can  be  generalized  to  the  more  complicated 
systems.   Also,  we  will  learn  what  restriction  must  be 
placed  on  the  coefficients  a^   and   p^   (v  =  0,1, ...,k) 
in  order  to  insure  convergence  of  the  difference  scheme 
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k  k 


(2.3)         5    a.Y  ^.   =   -L  h^"^  /    p.^ 

^   ^  ^^ — ;^   1  n+1         O       4 ^   1 


-M  -  -■-     -    1=0  -^^^-"^ 


where  we  recall  that  \x   =  h,      0<7<1. 

The  method  of  finding  a  solution  of  (2.3)  will  be 
that  of  Taylor's  series  or  small  perturbations.   We  can 
prove  the  following  theorem: 

THEOREM  2.1.   The  solution  to  the  difference  equation 

k  k 

4;^  1  n+1       o     4^  ^1  n+i 


where  L   is  a  positive  Constant  greater  tnan  zero  is 


(2-^)        \     '     Co«o^=l^l-  •••  +\-ld 


where  C  ,C-,,...,C,  .,   are  constants  and 
o   1      k-1 


(2.5)  Cj  =  Cj„  + 


L  p""'(c,.o)       ° 


l/m      2-27 


m 


+   Oih    "'     ) 


C.   is  a  root  of  p(C)   of  multiplicity  m  >  1. 
^jo  '^  — 

Proof: 

Case  1.   m  =  1.   We  look  for  a  solution  in  the  form 

of  a  power  series  in  h" 


1-7 


Let  Y  =  (C.   +  qh  "-^  +  •  •  •  )^  where  we  desire  to 
n    ^  ^jo   ^ 

find  the  value  of   q.   Substituting  this  expression  into 
(2.3)  gives 
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k 
1=0 


k 

y~  (a.  +  h^-^L  P.)(C.   +  qhl->'  +  ...  )^"'^ 
i^   1        °^    JO 


Dividing  both  sides  of  the  equation  by  Y^  and  using 
the  binomial  theorem,  we  find 

k 


The  first  term  on  the  left  is  zero  because  p[C^    )    =  0, 

l-y 
If  we  now  equate  coefficients  of  h   ,      we  obtain 


k 


y        (L  P.C^   +  ia.C^"^q) 


0 


or 


^  .  Vifjol   ana   p"(C.  )   /  0 

because  m  =  1. 

Case  2.   m  >  1.   Here   plCj^)  =  0   and   P^^ C^jq)  =  0. 
i  =  1,2, .. . ,m-l. 

In  this  case  we  look  for  a  power  series  expansion  of 
the  form 
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C.^  +  qh(l->')/"^  +  (7(h(2-2^^/"^) 


(a^.h^-Vi)(C^o-^^^;V'-^^^- 


i(i-l)  .i-2  2  2(l-7)/m  +  ...  + 


l(i-l)(l-2)---(i-(m-l))  .i-m  m  1-7  +  ...  )   ^  q. 
ml  jo  ^ 


Collecting  terms,  we  have 


(m-l),/.   N 
_    ,  P     ^'-jo^  ^m-l^(m-l)(l-7)/ni  ^ 

(m-l)l 

(m) / ^      \ 

[6(C.  )L  +- IJO-  q^Jh^-^  +  ..•   =   0. 

^  JO   o      m!      ^   ■' 


All  terms  but  the  one  in  brackets  and  those  following  it 

vanish  because   C  •   has  multiplicity  m.   Hence,  equating 

1-y 
the  coefficient  of  h  ^   on  both  sides  of  the  last  equa- 
tion gives 


°"/":.io'i-o 
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or 


That  is,   m  distinct  values  of  q  are  given  by  the  last 

l-y 
formula.   We  therefore  see  that  if  h  ^   is  sufficiently 

small  the  roots  C,  .   are  all  distinct  and  the  theorem 

J 

follows.   Q.E.D. 

A  simple  example  will  suffice  to  illustrate  the  fact 
that  some  well  known  difference  methods  must  be  excluded 
if  we  are  to  guarantee  convergence. 

Example  1.   Consider  the  equation  M-y'  =  -y  with 
initial  conditions  y(0)  =1.   We  choose  Simpson's  rule 
formula  as  our  difference  method: 


p(0  =  C   -  1;   the  roots  are   C  =  ±1«   Using  equation  (2-5) 


C.        =       1   -   h^'^  +  ... 
o 


and 


C,  =  -l-ih^-^. 


The  general  solution  is 
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Y^   =  C  (  1  -  h^-^   +  ...)^  +  C^(-l  -  ^  hl->'  +  ...)^. 
no  1      p 


However, 

»n     ,  -I  nH  , ,  ,  1  ,3.-7  ,     \n 

as  i.  -+  0  and  x  =  nh  by  (1.11 ).   Therefore,  since  we 
cannot  insure  that  C,  =  0   (due  to  rounding  errors  -  see 

later  discussion)  Simpson's  rule  is  divergent  for  boundary 

y 
layer  problems,  when  ix  =  h  . 

It  is  quite  clear  that  roots  of   p(0   which  lie  in- 
side the  unit  circle  will  not  cause  any  problems  with  re- 
gard to  boundedness  of  the  solution  of  the  difference  equa- 
tions.  However,  the  simple  roots  on  the  unit  circle  may 
lead  to  divergent  methods,  as  the  previous  example  illus- 
trates for  the  case   C  =-1.   C  =  1   is  always  an  acceptable 
root  because  by  the  consistency  condition   (p'(l)  =6(1)) 
and  by  equation  (2.5)  we  have  that 


C.   =   1  -  L y-^  +  (5(h2-27). 


o  o 


We  can  easily  see  that  the  condition. 


(2.7)  ^-^^^   <   0, 

p'(-l) 

ensures  that  ^  -*■  0     as  h  -*■  0  for  C-,   =  -1' 
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For  the  remaining  roots  on  the  unit   circle  we  proceed 

i  Q 
in  the  following  way:   Let  one  of  the  roots  be  Z-      -  e 

Of  course,  another  root  must  be  given  by  e~   .   Since  we 

have  assumed  that   d(0   and  p(C)   have  no  common  factors, 

6(e^^)  ^   0.   Let 


Pi  +  i^i 


d(ei^) 


p  (e   ) 
where  p.,   and  q,   are  real  numbers.   We  write 

C   =  (e   +  L^(Pi  +  iq;^)^  ^  +  ...)  . 

Taking  absolute  values 


,1-7 


l^^l    =    1^1^     =     [cos^e   +  sin  e   +  2Lq(p^cos   6   +  q^sin  0)h- 

.2,       2   ^        2s,  2-27   ^     <v^u2-27,Tn/2 

=     [1   +  2L^(p^cos0   +  q^sin   0)h^"^   +    (?(h^"^^)  ]  ^/^ 

Therefore,  a  sufficient  condition  for  convergence  in  this 
case  is 

(2.8)'  p,cos  9   +  q  sin  0   <   0. 
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We  can  state  the  following  useful  definition: 

DEFINITION.   A  difference  scheme  characterized  by 
the  polynomials   p{C)   and   6(C)   will  be  called  M.-stable 
if  in  addition  to  the  conditions  of  ordinary  stability, 
the  following  conditions  are  satisfied  for  those  roots 
of  p(C)   that  lie  on  the  unit  circle: 

1.  If  -1  is  a  root  we  require  that 

(2.7)  ^■^-=^   <   0. 

p'(-l) 

2.  Let  the  remaining  roots  on  the  unit  circle, 
excluding  C  =  1,   be  denoted  by  e  -'-,  e  2^  . .  .  ^  e 
Let 

ie . 

p.  +  iq.   =   -  -^^^-re^      J  =1,2,...,^. 

p'(e   J) 

We  require  that 


(2.8)     p  .cos  e.    +   q.sin  6.      <      0  j  =  1,2,...,^ 

J     J     J     J 


We  immediately  notice  that  conditions  (2.8)  really 
only  contain  i/2     conditions  because  we  are  dealing  with 
complex  conjugates.   That  is,  if  e  ^  is  a  root  so  is 

-ie 
e    .  But 
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and 


1  ■  -^1        ,   10, 
p'(e  h 


10^  -10-, 

„   +  ia    -     6(e   ^)   _     6(e    h 

p'(e   2)       p.(e    1) 


Pi  -  iQi 


where  we  have  let   0p  =  -0, .   Hence, 


Pi   =  Pg      and     q^   =   -qg. 


Now, 


PgCOS  ©2  +  qgSln  0^   <   0 


Is  equivalent  to 


p  cos(-0,)  -  q-,sln(-0,) 


<   0 


or 


p,cos  0^  +  q-,sln  0,   <   0. 


Therefore,  the  two  conditions  are  equivalent. 
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The  conditions  of  [x-stability  can  be  thought  of  as 
merely  conditions  on  the  coefficients,   3  ,   v  =  0,1, ...,k. 
That  is,  given  a  polynomial  p(C)   that  satisfies  ordinary 
stability,  we  attempt  to  find  a  related  d(C)   that  satis- 
fies the  consistency  conditions  and  the  ix-stability  con- 
ditions.  A  few  such  examples  will  clarify  this  point. 

o 
Example  2.   Let  p(0  =  C   -  1;   the  roots  are   C  =  ±1^ 

The  general  form  of  6(0   is   6(C)  =  ^^   +  ^^   +  P^* 

By  consistency 


(2.9)       p'(l)   =  2   =   6(1)   -  ^2  +  ^1  +  P 


By  condition  (2.7)  we  have 


or 


p'(-l)  -2 


o 


(2.10)  ^2  ■  ^1  "^  ^o   ^   °* 

Combining  (2.9)  and  (2.10)  gives 

Pl  -  Po'  -  Pi  +  Po  ' 


or 


(2.11)  ^1  ^   1- 
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Thus  the  inequality  (2.11)  is  equivalent  to  the  condition 

'2' 


2  .        2 

that  the  polynomials   p(C)  =  C   -  1   and   d(0    "^    '" 


P  C  +  P   satisfy  the  condition  of  M.-stability. 

We  now  see  why  Simpson's  rule  was  divergent.   In 


Simpson's  rule   P-,  =  ^,      which  does  not  satisfy  (2.11). 

Example  3 .   Let   p(C)  ==  C   -  1.   The  roots  are 

C  =  ±1  and      C  =  ±i.   6(C)  =  ^i^^'^    +  ^^^    +  ^2.^   +  ^iC  +  ^q' 
By  consistency, 

(2.11)*    4   -  p'(l)   =  6(1)   =  p^  +  p  +  Pg  +  P3_  +  3^. 


Condition   (2.7)    implies 


(2.12)        -^i^il     =      -  ^   O,    -   P      +  P      -    3      +  3    )      <      0, 
p'(-l)  443210 

Combining  (2.11)  and  (2.12)  gives 


(2.13)  2   >   P^  +  P^. 

For  the  root  t,   =  \.     we  have 

P   (1) 

or 

^1   =  -\^^\-   ^2  ^^o^     Pi   ^   -T  ^^3  -  ^1^- 
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We  require  by   (2.8)    that 


or 


q-^      <      0 


Pif    -    ^2    +  ^      "      0- 


Combining  the  last  expression  with  (2.11)  yields 


(2.1^^)  h      >     p^  +  2^2  +  P^. 

Equations  (2.1^)  and  (2.1^-)  express   i-L-stability  for  this 
example.  An  example  of  a  convergent  method  v;ould  be 


p(0   =   C^  -  1 


6(0   =   C^  +i  C^  +  C^  +i  C  +  1. 


2 


Check:   (2.1:>)  is  satisfied  because 


1    1 


2   >   2+2' 


(2.l4)  is  satisfied  because 


i|   >  I  +  2(1)  +|, 
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It  will  be  recalled  that  in  the  ordinary  differential 
equation  theory  stability  of  the  consistent  difference 
scheme,  i.e.  bo\indedness  of  its  solution  was  a  necessary 
and  sufficient  condition  for  convergence.   The  case  with 
boundary  layer  equations  is  a  bit  different;  that  is,   M-- 
stability  isn't  in  general  a  necessary  condition  for  con- 
vergence.  This  follows  because  we  have  not  considered 
the  cases  when 


(2.15)         p  .cos  0.   +  q.sin  9.    =   0. 
^  J     J      J     J 


Here  the  precise  determination  of  a  necessary  and  suffi- 
cient condition  seems  much  more  difficult  and  we  leave 
this  as  an  open  problem. 

In  order  to  be  complete  we  should  give  a  simple  exam- 
ple where 


p  .cos  0.      +     q.sin  9.        -        0 

16 
but  6{e        )    ^  0,      since  we  have  assumed  that   (3(C)   and 

p(C)   have  no  common  factors. 

Example  k.  Consider  p(C)=C^-C^+C-l.  The 
roots  are  C  =  1,  ±i.  p'(l)  =  2.  Therefore,  2  -  P^  + 
P   +  P   +  P   and  the  ii-stability  condition  implies 


28 


1   >   P^  +  ^2 


Suppose  instead,  we  let 


^1  ^  ^2   =   ^• 


Let   P^  =  lA   and   p^  =  3 A •   Also,  let   p^  =  p   =  l/2, 
Then   6(C)  =  -i  C^  +  |-  C^  +  ^  C  +  I  and 

6(1)   =2/0 

6(i)  =  -  ^  -  -^±7^  0,   p'(i)=-2-2i;   P  =  -  -^^   q  =  0 

6(-i)  =  -  ^  +  ;^i7^  0,   p'(-i)  = -2  +2i;   p  =  -  ^,   q  =  0. 


Not  only  can  examples  be  found  which  allow  p  .cos  9  .  + 

J     J 

q.sin  Q.      to  be  equal  to  zero,  but  we  can  also  give  examples 
J     J 

in  which  p  .cos  Q.    +  q.sin  9.   <   l/m  where  m  is  arbi- 

trarily  large.   In  practice  such  a  scheme  should  be  avoided 

because  the  convergence  of  the  scheme  would  require  a  very 
small  value  of  h. 

Example  5.   Let   p(C)  =  C^  -  1.   Then   6(1)  =  pp  +  ^   + 
P   =2  and  the  M--stability  condition  implies  P-,  <  1.   We 
choose 


6(-l)   =  ^2  -  P^  +  P^   =  l/m 
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1(1)  +  6{-l)   =   2^2  +  2P^   =   2  +  i. 


Let   P^  =1.   Then  p^  =^  and   ^1  =  1-^. 


Therefore 


C,   =  -1  -4:iiIhl->'  +  (^(h^^l-^^) 
^        p'(-i) 

=    -1  +4^1-^  +  Oih^^^-''^) 


is  a  root  of   p(C)  +  h"^  ^e(C)  -  0. 

Because  of  this  last  example  it  seems  wise  also  to 
consider  a  stronger  kind  of  stability. 

DEFINITION.   A  differnce  scheme  characterized  by  the 
polynomials   p(C)   and   d(C)   will  be  called  relatively 
stable  if  the  roots  of  p(C)  +  h  "■''6(C)  =  0  have  the 
property  that 

iql  1  ICJ  =  1  -  h^-^  +  (^(h^U-r))   i  =  i,2,...,k-i. 

For  a  relatively  stable  scheme  we  find  that 

6(-l)  <      _i 
p'(-l) 

and 

p  .cos  9.    +   q.sin  6.      <      -77. 
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These  conditions  are  more  stringent  thaxi  the  inequali- 
ties (2.7)  and  (2.8)  which  are  required  for  ix-stability. 

The  proof  that  ix-stability  implies  convergence  for 
the  linear  difference  scheme  (2.3)  requires  a  few  lemmas. 
Instead  of  using  subscripts  for  each  new  constant  that  is 
introduced,  we  will  use  the  notation  that  C   equals  a 
constant.   When  we  have  inequalities  like 


(C^  +  5)  <     C 


we  will  mean  that  a  new  constant  C-,   exists  such  that 


C^  +  5  ^  C  . 


Of  course,  the  fact  that  these  constants  are  independent 

of  h  for   0  <  h  <  h   is  to  be  demonstrated  for  some  h  . 

—  o  o 

LEMMA  1.   Let  the  consistent  difference  equation 

k  k_ 

>    a.Y  ,.   =   -L  h-^"^  /   P.Y  ,. 
4-^  1  n+i       o     4-3^  1  n+i 

be  M.-stable.   L   is  a  positive  constant  with  the  property 

that   -X  <  -L  <  -L  <  0.   Let  '^.   be  the  solution  of  this 
—    o  —  -*-! 

difference  equation  with  initial  conditions  $„  =  $-1  =  •  •  • 

=  \.2    =0'  \-l   =   (^k  -^  ^''Vk)'"- 
Then,  for  all  n. 
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n-1 

(2.16)  H  11,1  1  -^ 

1=0   ^      h-^  ^ 


for  h   sufficiently  small  and  where   C,  a  constant,  may 

be  chosen  independent  of  L  ,   i.e.   C   is  a  uniform  bound 

for  all  L   such  that   -L  <  -L  <  -L  <  0. 
o  —   o  — 

Proof:   The  solution  to  the  difference  equation  is 


I.  =  c^(h)?i  +c^(h)C^+  •••  +  Vi(h)?ti 


where   C  ^  C-, ,  •  •  • ,  C  _-,   are  linearly  independent  solutions 

of  the  difference  equation  and  are  distinct  if  h  is 

small  enough.   C  •  was  given  earlier  as 

J 

^i    =   ^io  +  ^^-  -Tiir^^  (-Lhi-^))i/^ 

J         JO  D^    ( C    ) 

^     ^   JO 


-f  ©(h^^l-:^)/^) 


where  p(C-  )  =0  and  m  is  the  multiplicity  of  C-  • 
By  Cramer's  rule,  we  can  write 
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C,(h)   = 


1 
C- 


^k-l   ^k-1 


0 


1 


'1 


^j+1   ^k-1 
^j+1    ^k-1 


e^  ,K-i...^.-i  (.^..i-r^A)-^  cti-  = 


k-1 
'k-1 


k-1 

2 

k-1 


,k-l 
'k-1 


The  numerator  we  denote  by  D.,   the  denominator  by  ¥ 

J 

(the  Wronskian) .   ¥e  observe  that  D.   and   W  reduce  to 
Vandermonde ' s  determinant . 


w    =    TT(q  -  c.) 

3<±  ^ 


,Djl 


(a,  +  h 
k 


'■v.> 


-1 


t<S 
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Therefore, 


C.(h) 


'  ^  k        ok    ' 


k-1 

i=0 

iq- 

■^J 

^M 

A  power  of  h  is  the  leading  term  of  a  difference 
expression  C,.    -   C,  .      only  when  C.   =  C-  •   Let  us  assume 

^  1      J       "^  lO      JO 

that  C        =1  and   C-   has  multiplicity  m.. 
oo  ^jo  ^         J 

Thus 


'^J^^^'   -   Ul-7)/in.J(m.-l)  • 
h       J   "^ 

If  m.  =1   (^-p,   is  a  simple  root),  then  by  p.- 
stability  there  exists  a  constant  L  >  0  such  that 
ICJ^  1  1  -  Lh-*""^.   In  this  case 

n-1  n-1 

lC.(h)|  H  ICJ'  1  cH  (1  -  Lhl-^)i/2 
'J     i=0   ^         i=0 


-  '^  -  (1  -  Lhl-^)^/2  ] 


<   c 


-  h^-^ 


for  h  sufficiently  small.   If  m.  >  1,   then  by  M-- 

J 

stability  for  h  sufficiently  small   [C-l  ^  r  <  1.   Here 

J 


y-^ 


OjWlgujI^  <_      ix-.)l..,-x.M,  g-^^ 


^ 

C        1    - 

^1-7    1    - 

-   r 

< 

C 

Com'bining  these  results  gives  (2.l6).   Q.E.D. 

Remark:   If  we  allow  y     to  be  equal  to  zero  and 
L   to  be  a  positive  or  negative  constant,  then  (2.l6) 
becomes 

n-1 

C 


I?.  1   <  -^   0  <  X   <  1   and  y   =  0. 
i=0   1      h     —  n  —         ' 


The  proof  is  similar  to  that  of  Lemma  1  with  the  use  of 

±Mx 
lim  (1  ±  Mh)^  =  e   ^  1  C 
h  -^  0 
x^=nh 

instead  of  (1.12) . 

Using  the  same  conditions  as  in  Lemma  1,  we  now  prove 
LEMMA  2. 


\\\  t  m  ^ii-X-D/m  +  '1  -  ^^'"''"' 
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for  h  sufficiently  small,  where  m  equals  the  maximum 

multiplicity  of  the  roots   ^  .   and  r  =  7^(1  +  max  |C  •  I  )  "=  1' 

JO  2     Kj^hl-^^ 

L  and  ^  are  positive  constants  independent  of  the 

particular  L   chosen. 

Proof:   In  the  course  of  the  proof  of  Lemma  1,  it 
was  shown  that 


In     '     Vo+<=l?l  +  ---  +Vl^k-1 


where 


iCjl  1       (l-yl(m.-l)/m.   "   l?Jol  ^  1 
h       J     J 


,0^1  1     C  If   ICjJ  =1. 


Also,  by  ix-stability  for  h  sufficiently  small 


Cjl  i  r  If   UjJ  <   1 


Cjl  <     1-   Il.l-''   if  ICj.l  =  1 


where  L  >  0. 

From  these  inequalities  the  lemma  easily  follows. 

We  complete  this  chapter  by  proving  the  following 
theorem  about  the  difference  method  associated  with  the 
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equation  (2.1);  i.e.   M-y'   =  -Ly+m(x). 

THEOREM  2.2.   Let  the  consistent  finite  difference 
equation 

k  k 

(2.17)   /   a.Y  _L.   =  h^"^  7   p.(-LY^.  +in(x^.))  +R 
^    '    4-^  1  n+i  4-^  ^1^   o  n+i     ^  n+i      n 

where  R   is  the  round-off  error,  satisfy  the  following 
conditions: 

a)  R^  =   0(h^"^^)   as  h  ^  0; 

b)  The  finite  difference  equation  is   M--sta.ble; 

c)  in'(x)   is  continuous  for  0  ^  ^  Jl  1^" 

d)  lY^  -  y(x^)l  <  Th-^"^   for   i  -  0,l,...,k-l. 

Then  the  difference  equation  (2.17)  is  convergent,  i.e. 
|Y  -  y(x  )  I  <  Ch-'-"'*'  for  n  =  0,1,2,  .. .   and  0  <  x^  <  1; 
C  is  independent  of  h. 

Proof:   The  proof  is  divided  into  two  parts.   We  first 

prove   Sj^  =  Y^  -  y(^n^  "  (7(h"^'^)   for   0  <  n  <  N^  - 

l-y 
—^ where   r  was  defined  in  Lemma  2  as  r  =  -p  ( 1  + 

max|C-  !)•   In  the  second  part  we  prove   ©^  =  0{h  ~ ^ )      for 

N  <  n  <  N  where  x,,  =  1. 
O  —    —  N 

First  we  note  that 


^iy"   =  -Ly'+ra'(x)   =  -L  (-Iioy_LJIl(2Sl)  +  ni'(x) 

O  O         ^l* 


or 


57 


h^y" 


,2(1-7), t2      r        f      \\      ,     u2-7  I,   V 

h  ^   ^  (L  y  -  L  m(x))  +  h   '^m  (x) 


=   (^(h^d-/)) 


where  we  have  used  statement  c)  of  the  hypothesis. 

The  exact  solution  y(x)   of  (2.1)  will  now  satisfy 
the  difference  equation 


k 


(2.18) 


i=0 


a. V  ,  . 


h 


1-7 


i=0 


p. (-L  y  ^.  +  m(x  ^. ) )  +  T 
1^   o^'n+i     ^  n+i      n 


where  y^  =  y(x^)   and  T  =  Cf{h^^'^  '^'). 

Subtracting  (2.l8)  from  (2.17)  and  letting  e. 


^4    -  y-  f      we  obtain 


(2.19)  IIc,  e^^.   =  -h^-'-L„IIPie„^.  +  (y(h^(l-''>) 


1^ 

i^ 


Let  the  right  hand  side  of  (2.19)  be  denoted  by 


n 


The  solution  to  this  difference  equation  is  given  in 
Hull  and  Luxemburg  [3]  as 


JZ   ^n-i-l  ^i 


(2.20)     e 


n 


i=0 

6 
n 


+     0  n  >  k 

n      — 


n  <  K 
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where  g   is  defined  by 


k_ 

=   0 


ll^i 


i=o  -  ^^-^^ 

and  g^  =  g^  =  .  •  •  =  g^_2  =  0,  g^_-^   =  a^'-^;    and  where 

k-1   k-i-1 

0        =     /        (  /   a,   .g„  ,,..-,  )e.   for  n  >  k  . 
n     4rQ  ^  ^     k-j^+k-i-j-1  i        - 

Consequently,  setting  maxlg^l  =  g  <  °°,   hy  M--stability 

k-1 


0  I   <   kag  X^  le.l   <   k^agTh^"^   <   K-^h' 


by  condition  d) . 


o 
Here  a  =  max  [a.  1   and  K-^  =  k  agT. 

Letting  P  =  max  \^.\      we  obtain  from  (2.19)  and  (2.20) 

n-1 
i=0 


le^l   <  h^-^PgL^le^l  +  h^-^pgL^(k+l)  2_  l^J 


n 


(h2(l-^))  +  10  1. 


+  SZ_  0 ,-n' 

i=k 

Now  if     h  _^  h        where      h^'-^pgL^  <    1      and      0  <   n  j^  N^   = 

^?  ^^'\        then 
Jin  r      * 

n-1 


(2.21)  le^l      1     h^~^  A  y~   le.l    +  K^h^ 
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where 


PgL  (k+1) 
A  =  2__ and  K, 


l-h^"^pgL 


o 


where  Kg  is  a  bound  for  R^  and  T^. 

From  (2.21)  and  a  simple  induction  (see  Henrici  [2] 
p.  2kk)    one  can  show 


le^l   <   K^h^-^(1  +  Ah^-^')^ 


(2.22) 


^  ,1-7  nh  ^A 
<  K-,h  '^  e 
—   P 


^  ,1-7  A(h^~^  In  h^~'^)/Jin   r 
<^  K,h  ^  e  ^ 


<   ^(hl-^) 


for   0  1  n  <  N^ 


For  N  <  n  <  N  we  consider  the  difference  equation 
o  —   — 


i=0 


^-1   -^  h'"Vi^^n+i  =  ^n  -  ^n 


with  the  solution  in  the  form 


n-k 


>   f  .  ,(R.  -  T. )  +  ^     n  >  k 
4-^  -^-1-1^  1    1     n      — 


(2.25)    e^  =  < 


^- 


n 


n  <  k 


^0 


where  ^   satisfies  the  difference  equation 

k 
>    (a.  +  h-'-"^L  p.  )^ 


1=0   -       -  -  ^^+^ 


with  initial  conditions 

k-1  k-i-1 

Y   =  7~~  {  J~   («!   •  +  h^"^p,   -L  ))?  ^1  •  ■   -,)e. 
n     4-^   -^^  ^  k-j       ^k-j  o'  ^n+k-i-j-1  i 

for  n  >  k. 


Since  here  we  are  interested  in  the  case   n  >_  N  , 
we  can  write,  using  Lemma  2  and  condition  d) 

k-1 

N. 

h" 

=   (^(h^"'*')  for  N  ^  n 


<   k^Ca  +  h^-^3L^)f(^^  +  (1  -  Lh^"^)'°)Th^"^ 


o  — 


N     N  in  r     „      ,1-7  . 

,          CO        ^n  h   '^  -ul-y 

because  r=e       =e  =  h        . 

By  Lemma  1  it  follows  that 


kl 


n-k  n-k 

TZ   I5n-i-ll|l'i    -  Til      1     H   lln.i.ilC^lh^'l-''') 

i=0        ^i   ^   -^        ^  1  i^O        "   1    -L 


Combining   these   results  with   (2.23)    gives 


(2.24)  e^      =      a{h^~^)      for     i\in<N, 


Q.E.D. 
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III.  OPTIMAL  METHODS 


The  criterion  of  ix-stability  leads  to  seme  interes- 
ting results  with  regard  to  optimal  methods.   By  optimal 
methods  we  mean  given  a  polynomial  p(C)   of  degree  k 

and  satisfying  p(l)  =0  what  is  the  maxim\im  order  that 

1  — "y 
can  he  achieved  in  the  operator  L[y(x);h  ~^]      by  choosing 

a  suitable  polynomial   6(C) •   For  the  ODE  case  the  order 

of  a  stable  operator  when  k  is  odd  cannot  exceed  k  +  1 

and  when  k  is  even  cannot  exceed  k  +  2.   We  now  briefly 

describe  the  theory  of  optimal  methods  as  given  by  Henrici 

[2]. 

We  associate  with  the  polynomials  p(C)   and   6(C) 

the  function 

(3.1)  MO   =  -^  -  6(C) 

where   ^n  C   is  made  single-valued  by  cutting  the  complex 
plane  along  the  negative  real  axis  and  setting  in  1  =  0. 
Using  (3'l)  and  a  few  simple  properties  of  analytic  func- 
tions, we  can  prove  the  following: 

THEOREM  3.1.  The  difference  operator  L[y(x);h] 
associated  with  the  polynomials  p(C)  and  6(C)  has 
the  exact  order  p  if  and  only  if  the  function  (t)(C) 


^3 


has  a  zero  of  the  exact  order  p  at   C  =  !• 

Proof:   See  Henrici  [2]. 

THEOREM  3.2.   Let   p(C)   be  a  polynomial  of  degree 
k  satisfying  p(l)  =0  and  let  k'   be  an  integer, 
0  _<  k'  j<  k.   Then  there  exists  a  unique  polynomial   6(C) 
of  degree  <     k'   such  that  the  order  of  the  associated 
difference  operator  is  at  least  k'+l. 

Proof:   See  Henrici  [2]. 

m     is  analytic  at  C  =  1,   so  suppose  we  let 


InX 


fii|  =  0^+C^(?-l)  +C2(C-1)2  + 


We  choose 


6(0   =  C^  +  C^(C-l)  +  •••  +  c^,(C-i) 


k 


By  the  theorems  we  have  that 


HO     =    Cp^^(C-i)P  +  ^((C-i)P^^) 


An  example  will  illustrate  the  technique  described  : 

Example  1.   Let   p(C)  =  (C-l)(C-^)   where   A   is 
real  aind  -1  <  A  <  1. 


4if 


P(C)  _  (C-l)[i-A  +  (C-i)] 
ITTT  "     iinLl  +  (4-l)J 


1-A  +s^{/:-i)  +^(c-i)2  .(^\c-i)^  +(5-((c-i)'^) 


Therefore, 


6(C)  =  1-A  +5zA(C-i)  +  i|^  (c-i)2 


1+5A    2-2a  .   3+A  >.2 
12      3   ^    12  ^ 


Note  that  for  7\  7^  -1   the  order  is  3   and  for  A  =  -1 
the  order  is  4.   Thus,  in  the  standard  ODE  case  it  is 
sometimes  possible  to  achieve   p  ==  k  +  2.   We  investigate 
this  further. 

We  introduce  a  new  variable 


(3.2)  >=   =  f^   °r   ?   =  T?| 


eind   two  new  polynomials   of   degree      <      k 


(30)  r(z)      =      (1^)^   p(i±|) 


(3.^)  s(z)       =      (^)^   6(^) 
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Henrici  [2]  shows  that  {J>.J>)    can  be  written 


2-'  +  • • •  +  -k-l^""'  +  -k^" 


(3 -5)    r(z)   =  a.,  z  +  a^z^  +  •  •  •  +  a,.  ^  z"  ""  +  a,.z 


where  the  coefficients  are  real  and 


(3.6)  a^  >  0 


(3.7)  a^  1  0      M-  =  2,3,  ...,k. 


We  now  consider  the  fimction 

(3.8)    p(z)   =  (1^)^  <i>(l±|)   =  -iM-  -   s(z). 

l-z 

The  function  p(z)   has  a  zero  of  order  p  at   z  =  0 
if  and  only  if  <|)(C)   has  a  zero  of  order  p  at   C  =  1 
and  thus,  by  Theorem  3'1,  if  and  only  if  the  operator 
associated  with  p(C)   and   6(C)   has  order  p.   Conse- 
quently, if  the  operator  has  order  p,   then 


(3.9)  s(z)      =     b^  +  b^z  +  bgZ^  +  ...   +  bp_^zP""^ 


where 


<'-^°>      itrgi?!  =  b„.v-v^- 
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The  degree  of   s(z)   cannot  exceed   k;   therefore,  a 
stable  operator  with  p  >  k  +  1  exists  only  if  ^i,  ,-1 
•••  =  b   -  =  0.   We  use  the  following  property  of  the 
in  C 

(3.11)        — ^  =  C^  ^C^z^  +c^z^  +  ... 
^^'  1^ 

where 


(3.12)  Cg^   <   0,     V  =  1,2, . 


If  k  is  odd,  then 

(3.13)      \+i   =  V^  ■-  C^a^^_2  4-  ...  ■+  C^^^a^. 

Using  (3.6),  (3.7)  and  (3.12),  it  follows  that  h^^^  <  0. 
Hence, 

THEOREM  3.3.   If  '^     is  odd,  the  order  of  a  stable 
operator  cannot  exceed  k  +  1. 

If  on  the  other  hand   k  is  even,  then 

(3.1^)  Vi    ^   '^2\^h\-2  -"  '"  ^  V2  • 

A  necessary  and  sufficient  condition  for  ^v+n  =0  is 
that  ap  =  a.,,    =  •••  =  a,  =  0.   That  is,   r(-z)  =  -r(z). 


^7 


By  stability  the  real  parts  of  the  roots  of  r(z)   must 
be   <  0.   Since   r(z)   is  an  odd  f-unction  here^   r(z) 
ceinnot  have  roots  with  negative  real  part.   We,  therefore, 
have  by  (5.2)  that  all  roots  of   p(C)   lie  on   |C|  =  1. 
Since  a,  =  0,   the  degree  of  r(z)   is  k  -  1.   Thus   -1 
is  a  root  of  p(C)   (otherwise  r(z)   would  be  of  degree 
k).   Using  (3.6),  (3-7)  and  (3.12)  again 

(3.15)    \+2   =  S^k-1  +  '^6^k-3  '"•••"'  ^k+2^1 

is  negative,  we  find  that  the  order  of  a  stable  ODE  opera- 
tor cannot  exceed  k  +  2. 

THEOREM  3.^.   A  necessary  and  sufficient  condition 
for  p  =  k  +  2  is  that  k  be  even,  that  all  roots  of 
p(C)   have  modulus  1,  and  that   6(C)   be  determined  by 
(3.5). 

We  are  now  prepared  to  im.pose  the  condition  of  M-- 
stability  and  investigate  the  changes  in  the  above  theory. 
It  is  quite  clear  that  if   C^^  =   1  and   K-^l  "^  1 
(O  =  1,2,  . .  .,k-l),   then  we  can  use  formula  (3-l)  "to  obtain 
a  difference  operator  of  order  p  =  k  +  1.   We  also  have 
the  following  result: 

THEOREM  3.5.   All  stable  methods  of  order   p  =  k  +  2 
are  not  (x-stable. 
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Proof:   Since   C  =  -1  must  be  a  root  of  p(C)   in 


order  that  p  =  k  +  2,   it  will  suffice  to  show  that 

_d(-l) 

P 

form 


■^1 — YY  <   0  cannot  be  satisfied.   Since   p(C)   has  the 


10,      -10.  ^^        ^  "^ 


p(a   =   (C-l)(C+l)(C-e   ^)(C-e    l)...(C-e   2)(C-e    ^) 


-K-2       -V2 


p'(-l)  =  _2(-l-e   ^)(-l-e   l)...(-l-e   ^  )(-l-e    ^  ) 


<   0. 


Therefore,  we  need  only  show  that   6(-l)  <  0.   Since  by 
(5.2)   C  =  -1  corresponds  to  z  =  00,   we  have,  using 


d(-l)   =  2^  lim-^^ 
Z-+00  ( l-z) 

,       b  +b-,z+«  •  •+b,  z 
=  2^  lim  -^-^^ ^-^^ 

-   Z— >-oo        (1-z) 

2"  >^'  \ 


(-1)"  k! 


1a- 

=  2  b,     by  L' Hospital's  rule 

K. 


Since  k  is  even. 


^9 


\      =  ^2^k-l  +  %^).-J>    +  •••  +Ck^l  • 


Using   a^  >  0,  a^  >  0   (v  =  2,3,-..),   and  C^^  <  0, 
we  obtain  b,  <  0  and  therefore   6(-l)  <   0. 

Hence,    i/"-,  \    >   0  and  the  result  follows. 

Immediate  consequences  of  this  theorem  are  the  following 
corollaries: 

COROLLARY  1.   If  k  is  even  and   p(-l)  =0,   then 
the  unique  stable  scheme  of  order  p  =  k  +  1  is  not 
M.-stable. 

Proof:   Since  k  is  even  and  C,   -  -1     is  a  root, 
there  exist  an  even  number  of  roots  of  p(C)   which  are 
real.   Denote  them  by  1, -1, A  , Ap, . . . ,A^  .   As  in  the 
theorem,  we  have 

p'{-l)   =  -2{-l-A,)  •••  (-l-^ps^   (positive  number). 

By  virtue  of  stability   |a.1  <  1.   Thus  p'(-l)  '^   0.   As 
in  the  proof  of  the  theorem  (3(-l)  <   0.   Thus,  the  scheme 
is  not  M--stable. 

COROLLARY  2.   If  k  is  odd  and  the  roots  of  p(C) 
are  A,  =1,  A  =  -1,  A_,...,Aj(-  where   |A.1  <  1  for 
i  =  J>,k ,  , . . ,]<i.      Then  the  unique  stable  scheme  of  order 
p  =  k  +  1  is  |J.-stable. 
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Proof:   By  hypothesis  it  easily  follows  that   p'(-l) 
>  0.   Since   k  is  odd 

\   =   Vk-1  '-h\-j>  +  •••  +  \-A  • 

Using 


Cgy   <   0    V  =  1,2,.. 


and 


^   ^ 


0     [I    =   1,2, 


we  obtain  b,  <  0.   But  b   =  0   is  equivalent  to   6(-l)  =  0 

K  —  K 

((3(_l)  =  2^3,  )   and  we  have  excluded  the  case  of  common 
factors  of  (5(C)   and  p(C).   Therefore,   d(-l)  <  0  and 


^^  ^^   <   0.     Q.E.D. 
P  (-1) 


In  concluding  we  say  that  the  best  precision  that  we 
can  hope  for  is   p  =  k  +  1.   From  Corollary  1  we  see  that 
this  value  of  p   is  not  always  obtainable.   In  general 
when  p(C)   has  roots  on  the  \mit  circle  and  k  is  odd, 
we  must  test  separately  whether  each  complex  root  on  the 
unit  circle  satisfies  the  conditions  of  la-stability. 
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Example  2.   Let   p(C)=C^-C^+C-l.   The  roots 
are   l,±i,   and  the  M--statiility  condition  reduces  to 
1  >  P   "^  Pp*  '^^^   unique  optimal  method  for  this  particular 
choice  of  p(C)   is  characterized  by 


d(C)   =  3^(5C^  +  7C^  +  7C  +  5) 


where 


^1  +^2   ^  i  +  i  =  i  ^'   1- 


Therefore,  this  optimal  method  is  not  M.-stable. 


IV.   BOUNDARY  LAYER  EQUATIONS  OF  THE  FIRST  ORDER 
IN  ONE  UNKNOWN 


In  this  chapter  we  discuss  the  numerical  treatment 
of  the  equation   (1.6)   M-y'  =  f(y,x)   over  the  interval 
[0,1]   with  initial  conditions  y(0)  =  a.   We  will  ass\ime 
a  slightly  more  general  condition  of  positive  stability 
than  the  one  given  in  Chapter  1.   In  addition  we  require 
that  the  real  parts  of  the  roots  of  the  characteristic 
equation 


Det  (||  (y(x,iJ.),z(x,M.),x)  -  Al)   =  0 


are  negative  and  boimded  away  from  zero  in  D  for  all  u- 
sufficiently  small.   For  the  equation  (1.6)  this  condition 
reduces  to  -r—  j<  -L  <  0  where  L  is  a  constant. 

We  will  naturally  assume  that  all  of  the  conditions 
of  the  hypothesis  of  Theorem  1.2  are  satisfied. 

THEOREM  k.l.      Let  the  consistent  finite  difference 
equation 


(4.1)        XI^.Y^+i   =  h^"^^   P.F.^.  +R 


i^  i  n+i  ^  i  n+i    n 
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where   Fv.+-4  =  •^^^n+i*"^n+' ^   ^^^  ^n  ^^   ^^^  round-off 
error  satisfying  the  following  conditions: 
a)   R^   =   (7(h2(l->')); 

h) -   The  finite  difference  equation  is  relatively  stable; 

^  -p   A  -p       ^  -p 

c)  -^f      -r—  and  — p—  are  continuous  and  bounded 

^    y       by 

( -E  <_  ^  <_  -L  <  0)   for   0  ^  X  :!  1   and   -oo  <  y  <  oo^ 

d)  |e^|  =  |y^  -  y. I  <  Th^"^  for  i  =  0,1, . . .,k-l  where 
T   is  a  positive  constant  independent  of  h. 

Then  there  exists  a  constant   C   for  a  given  value  of  h 
such  that   |e  I  <_  Ch  "■^  for  n  =  0,1,..., N  where 

Proof:   The  exact  solution  of  [ly^    ^  f(y,x)   satisfies 
the  difference  equation 

k  k 

(^.2)        y        a.y  ^.   =  h-*-"^  ^   p.f  ,.  +T 
1=(J  j_=U 


where  T   is  the  truncation  error  and   f  ,.  =  f(y(x  ,.), 
n  n+i    ^  -^  ^    n+i  ' 

X   .).   By  our  hypothesis  on  ^  and  y-  and  the  consis- 
tency condition  we  can  show  that  T   =  0(h  ~  '*')   as 
h-^0.   This  is  done  as  follows:   By  Taylor  series  we  have 


where 

X    <   0  ,  .   <   X  ,  . 
n  —   n+i  —   n+i 


5^ 


and 


where 


f 

-nil     =     y  ^.      =     y      +  ihy"(0    ^.  ) 


n     —       n+i     —       n+i 


Also, 


."    =   l_iflyjjil)    -    f    y'  +  ^x    _    V  .^  ^x 

dx^      ix  y    \x  \j.  2         m- 

M- 


Substituting   these   expressions    into    [k .2) ,    we    obtain 
k  k 


r       /        a .    +  hy       )        ( ia .    -   p . ) 


"      i-O      "  ^^     i^ 


^^  2  ■'^ 

+     y~"  a.ii^L      ..(0         )    _   ^2      Y~(3.(iy..(a'       )) 
4 — rT     1      2  n+1  4 — TT  '  1^    -^    ^    n+i 


=     T    . 
n 


Or   since     ^ a.    =0      and     /. {i-O..    -   p.  )    =  0     by   the 

i=0   ^  i=0    ^    ^ 

consistency  condition,  we  have 
n       4-rQ  2  j^27         ^-^^  n+i  '  n+i 

+  ^x^y^Vi^-Vi^  . 
h^ 
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=  Q(h"^"^^)    as   h^O. 


We  now  return  to  the  proof  of  the  theorem.   Subtracting 
equation  (4.2)  from  [k  .1)    and  letting   e.  =  Y.  -  y.,   we 
obtain 


i=0  -   '^'^  i=o 


K 


n    n 


h^'^H^-  If^^^  e  ^.  +  R   -  T 
4-3q   1  dy     n+i    n    n 


where  ||^+1  =||^^^+i  ""  Vi^^n+i"  ^n+i^ '^n-fi^ '  ^'^  ' ' 

-vrr  is  evaluated  at  a  point  between  y   .   and   Y 
oy  "^  "^n+i        n+i 

If  we  denote  the  right  side  of  (4. 3)  by  Q  ,  we  may  pro- 
ceed as  in  Theorem  2.2  using  the  following  bo\ind  for  O 


<   h-^"^  Ep  )~  |e  ^.  I  +  |R   -  T 
—  ^^ — -  '  n+i '    '  n    n 


i=0 
where 


P   =  max  IP^I   . 
i=0,l, . . .  ,k 
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It  will  be  recalled  that  in  Theorem  2.2  for   0  <_ 
1-y 

^  ""  N   =  -j^— where   r  =  o(  1  +  inax  I  C  .^  I  ) ,   it  was 

—   oi:.nr  d  \  tl       \  <1^ 

'  ^jo' 

possible  to  obtain  an  explicit  bound  for  e  .   If  we 

replace  the  L   in  Theorem  2.2  by  E,   we  obtain  the 

following  bound  for  e   in  the  non-linear  case 


e^l   <   K^h^-^(1  +  Ah^-^)^ 


1-7   nAh 


1-7 


<   ICh   ^  e 
(1^.5)  l^n'   -  S^    ®  for  0  <  n  <  N^ 


2  o       1      J    n      pg-L/lK  +  1; 
where  K_  =  — - — ^-- — = —   and   A  =  "^^   ^  — = 

^        1  -  h^"^PgE  1  -  h^'^PgL 

l_-y    = 

where   h   '^PgL  •<  1   and   h  <_  h  .   K   is  a  bound  for  R^ 

2 
and  T  ;   g   is  the  Green's  function  bound,  and  K^  =  k  agT, 

the  bound  due  to  the  initial  error,  where   a  =  max|a.|. 

We  now  find  an  interval  [0,Xp.  ]   where   -L,  <_ 

k-1,  k  (l-k  1   ) 

(if. 6)   h^-^XI  (L,  -  L.)  H  IPJIC  I     1^  I     <   1 

where  the   '   means  the  s\im  is  taken  over  those   i's 
corresponding  to  roots   IC-^I  =1  and  the  C . ' s  are 
defined  with  respect  to  the  difference  equation 
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k 
i:  (a.  +  hl-'-p.L^)!^^.   =  0 


with  initial  conditions  ^  =!li'^'*"'=!^vp=0   and 
$1,  -,  =  T — - — ^—   .   The  Green's  function  takes  the  form 

$„  =  c^c!^  +  c,C?  +  ...  +  c,  ,cP  n 
-HI      o^o    11         k-1  k-1 

and  by  relative  stability 

iql  1  ICqI    i  =  l,2,...,k-l 

where 

C^   =   1  -  L^h^-^    +  (y(h2(  !->')). 

By  a  close  examination  of  {h.6)    it  becomes  clear 

that  the  length  of  the  interval  [0,x,,  ]   depends  upon 

^1 
the  value  of 

k-1  ,      k 

Id  y~   le  . 


(^•7)  1  1  ZZ  |c.|  2_  le.|. 

i=0       j=0    ■' 


We  shall  see  later  how  to  choose  the  a.   and   S.   to  make 

1        1 

this  sum  close  to  1.   For  now  we  will  merely  assume  that 

the  sum  is  close  enough  to  one  so  tha.t  N  <  N-,  . 

1-  V~^     =      o    1 
If  we  now  add   h  "■^  /   .  „  p.L-,e  ,.  to  both  sides  of 

*■ 1=0      1   1  n+i 

equation   {^ -3) >    we   obtain 
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k 


(4.8)     JZ   («i  +  h^-^p.L^). 


e  ,  . 
•i=0   -        ^  ^  n+i 


1-7  y~  A  /f   ^  J+1- 


1^  1   dy  n+i    n    n 


i=0 


We  let  the  right  hand  side  of  (4.8)  he  denoted  oy     q  , 
and  we  note  that 


=      ,    Sf(y   .  +5  ,.e  ,.,  X  ,.) 
L-,  +   '"'n+i    n+i  n+i'   n+i 


^   ^ 


=    af(y^^.,x^^. )    ,2^ 
1   '^1  -^^  +^  Vi^n+i 


for  0  ^  n  <_  ]\^  -k 


<L,  -L,   +/^(e^.)      ,.    „,-^, 
—   1    1     ^^    n+i    ana   i  =  0,1,  ...,k 


^2 
where  we  have  used  condition  c),  i.e.,   — ^  is  continuous 

and  bounded. 

Now  it  follows  that 


k 

+  (^(h2(l-^))  4-{^(hl-^II  ,     ,2. 

j=0  '^i+j'  ^ 


0  1  i  1  ^1  -  ^ 


using  condition  a)  and  our  knowledge  about  the  behavior 
of  the  truncation  error. 
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The    solution   to  the   difference   equa.tion   (4.8)    is 


(4.10) 


n-k 


e        =     < 

n  ^ 


i=0 


n 


^      .    ,0.    +  ^  n   >   k 

-'-n-i-l^i  n  — 


n  <   k 


where  ^   satisfies  the  difference  equation 

_k_ 
> 

1^ 


II  («i  +  hl-^p.L^)^^ 


•n+i 


=  0 


with  initial  conditions  ^  =  l^-i  =  •  •  •  =  ^k-2  "  ^  ^''^'^ 

ik-1  -  (-k  +  h''Xh)"'- 

^   is  given  by 

k-1  k-i-1 


n 


i=0    j=0 


(a,   .  +  h   -^p,   .L-,  )¥  ,.  .  .  ,  )e. 
^  k-,1       ^k-,1  1  -J^n+k-i-j-1   1 


for  n  >  k 


However,  for  n  _>  N   we  have  by  Lemma  2  that 


N 


o 


n    N 
'■J-n'   —  -^  ^  ,  (1-7}  (m-lT/m   ^         '   '    o  — 


where  L,  r,  and  ^  are  uniform  bounds  for  a].l  L-,   between 
L  and  L,   and  where  m  equals  the  maximum  multiplicity 
of  C  .  •   Therefore,  using  condition  d) 


6o 


k-l 


l^n'  1   («  +  hl'^PL)X(lf^^^_J  +  ...  +  11^1)  YZ   lej 


i-0 
N 


o  .  ..  N 

h 


1  k2  (a  +  h^-^PL)©  V^  +  (1  -  Lh^-^)  °)Th^->' 


1-7 


(4.11)1^  I  <  K,,Th  '  for   n  >  N 


where 

N 

h 


Kj^   >   k2(a  +  h^-^PL)f(  -^  +  (1  -  Lh^-^)  °). 


By  using  Lemma  1  we  can  find  a  boimd  for  the  part  of 
(4.10)  that  is  due  to  the  round-off  and  tr-oncation  errors 

n-k 
(^•12)       II  IIn.i.il^(h2(l->'))   <   K^h--^ 

where   K_   can  he  chosen  as  a  uniform  hound  for  all  L 
5   _   _  -^ 

such  that   -L  j<  L  <  -L  <_  0. 

If  we  define  E   =  raax(  | e  I ,  1 e,  1 , .  .  . , | e  |  ) ,   then 
by  using  Lemma  1  again 

1=0 


A  bo\md  can  now  be  found  for  e   using  the  last  inequality, 
(4.9),  (4.11),  and  (4.12)  in  equation  (4.10): 
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ej  1  Z:ii;,.i.iikii  +  i\i 

k  n-k 

-L      j-^o        ^        ^1=0        nil 

H-  5  on 

k  n-k 

-L      j=0        '^      "   i=0  °        ° 


2 
6     n 


k-1,  k 


i=0  ^  ^      j=0        J        ^        l-lC^I 

n  4-  p  on 

where  the     uih^    -7'/^^      term  results   from  the   roots      It-    I 
<   1     and  Lemma  2,    i.e. 


\1A       <      1( ^- +    (1    -   Lh^-^)^) 

^       -  ^(l-7)(m-l)/m 

where      r   <   1. 

By  our  choice  of  L-,   and  L-,   (equation  (^.6))  and 
assuming  h  is  sufficiently  small  the  term  in  brackets 
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multiplying  E   will  be  less  than  a  <   1. 
Hence, 


(^.15)  le^l  1  aEn  "*"  ^5^"^"^  "*"  Ki^Th"^"^  +  KgE^  for  N^  ^  n  <  N^, 


If  E   =  le.l   for  j  <  N  ,   then  we  may  use  our  esti- 

n    '  J  '  '^  —  o'  "^ 

mate  (^.5)«   If  J  ^  N  ,  we  can  replace   |e  I   in  (^.13) 
"by  E   and  obtain 


(if.li^') 


E 


<   A  ^^-7   ^   "^ 


n  —  1-a 


K.,Thl->'_^KgE2 


1-a 


1-a 


for  0  j^  n  <_  N^ 


since  this  bound  is  larger  than  the  one  for   je  I   for 

0  <  n  <  N  ,   i.e.,  see  (4.5). 
—   —  o 

Geometrically,  we  have  the  following  picture: 


E   -  T— ^  E 
n   1-a  n 


K,,  T 


n 
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We  wish  to  show  that  E  <   d,   for  0  <_  n  <_  N^.   It  will 
then  follow  that 

"^1  -  1^  "^1  -  1?S    -   1-a    =  ° 


or 


-1  ./l  -  ^  ^—-1  Kgh^  ^ 
>/         (1-a) 


^1     


-2Kg_ 
1-a 


1  —  a 


The  proof  that  E   <_  d,   for   0  <_  n  <_  N-^   is  by  induc- 
tion.  Rewriting  the  difference  equation  (^O)  with  n 
replaced  by  n-k,   we  have 

For  h  <  h   where   |a  |  >  Ejp  jh  "'*'  we  have 

—   O  K.  K.    O 

1=0 


'"    -  Ki  -  ^Kl'-l'' 


^   K-k  -  Vki 


K\-^KK^-' 
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i  "V  \-l   +  Kg  h^'"-^' 


where 


^l_i ^^ iBI 1     ^ 


and 


lo^l-ElpJhl-/  -   8 


\-l  i  T"^ 


<  Th-^"'^  < 


d,   by  hypothesis.   By  induction  assume 


that  E  -1  ^  d,   and  h  is  so  small  that 


l^n'   -  ^  ^1  +  ^8  h^^^"^^   ^   ^1+^2 


and  therefore  E   <  d   +  dp.   But  now  by  (^.1^')  we  have 

that  E  <  d,  +  d„   implies  E  <  d.,  .   Therefore,  by  in- 
nl2^       n—  1  ^0- 

duct ion 

K_  h       K  Th 
(i^.lif)  E^   <  -^ -+4 ^  +  ^(h2(l-r))   for   0  <  n  <  N. 


We  may  now  repeat  the  argument  and  extend  the  interval 
to   [0,XjT  ]   where   N-,  <  N^.   We  find  an  interval   [x„  .v+t^ 

Xjj  ]   where  -L^  1  ^  jl  -Lg  and 
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k-1,  k  (1-K  I       ) 

(4.15)    IZ  hi-^(LVE,)  IZ  iPjl  icj  iql"^-^       ,\  I,.  I 


<   1 


where  the  C. 's  are  now  defined  with  respect  to  the 
Green's  function  associated  with  the  difference  equation 

k 

1=0 

with  initial^  conditions  %^  =  '^■^  =   '"    =  %^_2   "  °  ^"^ 
Ik-l  =  (°'k+^''''\^2)''- 


We  add  h^"^  )   , _^  P,  L^  e„^,   to  both  sides  of 


-k 

-i=0  ^i  ^2   "n+i 
equation  (4  O) : 


T~   (o^i  +  h^"^  Pi  L2)e^ 


i=0 


l.ry^..^^2-S"'^^n.i^^n-V 


h-'-W   P.( 


i^o  ^ 


We  treat  this  equation  exactly  as  (^.8)  was  treated 
hut  we  use  (if.l4)  for  the  initial  estimate.   The  accumula- 
tive error  for  the  first  two  intervals  is 


-r.    .^'^"^^I^I^'^"^^!^-"^!^*"^'''"'' 
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Ki|    K_   ,  _^  K),  2 


i  II  +  ir|]  ir|h^-'  +  (ttI)  Th^"''  -  ^(h^''-''') 


K,  2 
[1  -  (-^)  ]   K^   T_       K,  2 


^ir^irlh^"'  +  (ir|)  Th^-''  +(?(h2'^-^') 


0  <   n  <   No  • 


At  the  end  of  the  third  interval  the  error  will  be 

K,,  2 

Kp;    1_         Ki,   [1  -  (- 

E 


<  A  ,1-r  ,  A  '^  -  't^'  ^  A  ,1-r 


■"n  —  1-a  1-a       K^     1-a 

1  -irl 


X  —  3, 


li^C^.  (^)\h-  .^h^U-.)^ 


1    ^^ 


iTi: 


0  j^  n   <_  N^  . 


Finally  at  the  end  of  the  Ji      interval  we  obtain 

1— iril 


for   0  _f_  n  ;<_  N.  =  N.   Q.E.D. 
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Of  course,  the  bound  for  E   given  by  (4.l6)  may  be 
large  if  i      is  much  greater  than  1,  i.e.  if 

k-1,      k 

2 |C.  I  2 l^n-l   is  much  greater  than  1. 

i=0      j=0   J 

k 

We  therefore  wish  to  minimize  this  sum.   / |p.|   will 

J=0   J 

be  a  minimum  if  we  choose  all  p.  >_  0.   Since  the  roots 

J 

C.  ^  1     on  the  unit  circle  do  not  yield  ix-stable  difference 

schemes  of  higher  precision  than  those  roots  inside  the 

unit  circle/  we  will  exclude  such  roots  for  now. 

We  are  thus  left  with  the  problem  of  minimizing 
k 

|C  1  2 |P^-|.   The  following  corollaries  will  be  very 

°  j=0   J 

useful: 

COROLLARY  1.   If  the  only  essential  root  is   ^ .   =1, 

^  JO     ' 

and  if  a,  =  1  and  p.  >  0  for  j  =  0,1,  ...,k,   then 
J^-  J  — 

the  value  of  Z      in  Theorem  4.1  is  one. 

Proof:   Since  p .  >;^  0  by  consistency 

J 


k 

(1)"  =  ZI  IP. 


j=0   ^ 


=  (1  -  r^)(l  -  r^)    ...  (1  -  r^_^) 


where 


p(r^)   =0     i  =  1,2, . . .,k-l, 


68 


In  Lenuna  1  it  was   shown  that 


< 
o'      — 


l(V«l'(?o-^2'  •••  (V<k-1> 


1  +(:^(h^-y) 
p'(i) 


Therefore, 


'  I  7Z  IP 

O'   ^7— TT   '   , 


~   1 

j=o   "J 


and  i,      can  be  chosen  equal  to  1.   Q.E.D. 


COROLLARY  2.   If  the  only  essential  root  is   C-   =  1 


and 


(L  -  L)|C^| 


k 


°    '        'P  J   <  1 


L  +(7(h^"^)   j=0   J' 


where 


KJ   =   1  -  Lh^-^  +  (^(h2U-7)), 


then  the  value  of  i      in  Theorem  4.1  is  one. 

Proof:   The  proof  is  obvious  since  (4.6)  is  satisfied 

for  all  0  <  x^  <  1. 
—  n  — 
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We  will  consider  a  few  examples  to  illustrate  how  the 
value  of  i,     may  be  estimated  in  practice.   Consider  the 
differential  equation 

M-y'   =   -y{y  -  l){20x  +  lO);     y(0)   =  2. 


In  the  bo'undary  layer  for   all     0  <   h  _<_  h 


-30     1     f       1     -10. 

•J 


Example  1.   Suppose  Adam' s, method 

Yn+5  -  ^n-f2   =  ^  ^  ^F^^^  +  19F^^2  -  5F^^,  +  F^) 

is  used  to  calculate  the  solution  to  the  above  differential 
equation.   Here 


TL  IP,'  '  2i 


o  '      '■ '  ^1 '      24 

Since  f   is  monotonic  in  the  boundary  layer  L,  =  L^  =  L,  , 
Lp  =  L,  =  Lp,   etc.  ^^'     ~ — 1  ^  =  .8  where  we  have  let 
a  =  .8  <  1. 

L^  =  30  (1  -  .8  (|^)) 

=  30  (.^35)   =  13.1  ■ 
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Lg  =  30  (.^35)^  =  5.66 


L^  =  30  (.435)^. 


Since  Lp  =  5*66,  it  is  clear  that  ^  =  2  is  a  good 
estimate  for  the  value  of  Z      in  the  bo\indary  layer. 

Example  2.   Consider  all  methods  with  essential  root 
C  =  1  and 


C    <   2. 
o'   — 


We  choose  a  =  .8  again.   Then 


L^   =  30  (.6)^ 


L^   =  18 


L,   =  6.5  . 
3 


Thus  we  take  i  =  3   in  the  boundary  layer. 
Example  3.   The  same  as  Example  2,  but  let   |C^1  jf.  3  • 
Here 


L^  =  30  (.73^)^ 


^4  ~  8.7  • 
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Therefore,  we  choose  £   =  ^4 . 

The  above  calculations  were  made  inside  of  the 
boundary  layer;  outside  of  the  boundary  layer  we  make  use 
of  the  following  fact  from  the  asymptotic  theory: 


fy(y{x,tx),x)   =  fy(4)(x),x)  +  0{h'^) 


for 


-  §  h^  in  h^   <_  X  ^  1 

L 


In  our  examples   o(x)  =1   so 


f  (l,x)   =  -20x  -  10 

J' 


outside  the  boundary  layer.   Here  f  (l,x)   is  independent 
of  h  and  is  monotonically  decreasing  from  -10  to  -^O. 
Therefore,  in  Examples  1,  2,  and  3  we  must  increase  i     by 
2,  5,  and  4  respectively  in  order  to  include  the  entire 
interval   [0,1]  in  our  estimate. 

In  practice  it  was  observed  that  there  was  no  error 
build  up  outside  of  the  boundary  layer  region  for  ii-stable 
schemes.   Therefore,  the  above  estimates  are  to  be  considered 
as  absolute  maxima  for  the  value  of  Jl . 

A  few  remarks  about  Theorem  ^.1  are  in  order: 
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Remark  1:   The  same  proof  could  be  used  to  obtain 
bounds  for  tx-stable  methods  instead  of  the  less  general 
relatively  stable  methods,  but  these  M.-stable  methods 
would  require  a  much  larger  value  of  Jl . 

Remark  2:   Condition  to  (c)  could  be  weakened  by 

requiring  that  the  functions  ^,  -^     and  — ^  be  con- 

^  y  Sy 

tinuous  in  a  strip:   0  Jl  ^  :!  1  a-^d   |y  -  y(x)  |  <  t  where 
t  is  as  large  as  is  necessary  in  the  proof. 

Remark  3:   By  using  the  variational  equation  and  by. 
weakening  conditions  (a),  (c),  and  (d)  slightly  we  can 
gain  another  0^{h   ~^)      in  our  error  bounds.   This  is  done 
as  follows:   We  make  the  following  changes  in  the  hypothesis 
of  Theorem  ^ .1: 

(a-)   R^  =  ^(h^(l->'^) 

(c')   f(y,x)   has  two  continuous  and  bounded  derivatives 
for  0  j<  X  _<  1  and  -oo  <  y  <  oo. 

(dM   e^   =  ^(h^^^-^'^)   i  -  0,l,...,k-l 

Using  conditions  (a'),  (c')  and  (d')  we  derive  the  so-called 
variational  equation. 

Using  condition  (c')  and  (4.16),  we  may  write 

f(Y^,x^)  -  f{y(x^),x^)   =  fy{y(x^),x^)e^  +  aih^^^'^h  . 

Our  difference  equation  {^  .J>)    takes  the  form 
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k  k_ 


(^.17)   )   a.e  ^.   =  h-^"^[  )   p.f  (y(x  ^.),x  ^.)e  _^.  + 
^    '    A — ;r-  1  n+1         "■  4 — ?^  1  y^^^  n+i  '  n+i  n+i 


0[h^^^-'^h]    -   Cy"(x^)h2  +  ^(h^(l->'^) 


where  condition  (a')  was  used  for  the  round-off  error 

and  Cy"(x  )h   for  the  truncation  error. 
^  e 


h--^' 


If  we  let   e   =  ■■  ^ ^-  ,  equation  (^.17)  becomes 
n   ^1-7 


(1^.18)  IZaie^^i   =  hl'^'t  ZI  Vy(y(^n+i)'^n4-i^^n+i^ 


-  Cy"(x^)hl-*-^  +  a(h2(l-^h  . 

Since  h  ^y"{x  )   is  a  continuous  function  for  all  \i, 
we  can  write 


k 


h^^  Cy"(x^)   =  C  h2>'2_P-  y"(^n4-i)  +  ^(h) 
"  i=0   "^ 


Substituting  the  last  expression  into  (4.l8)  gives 

k  k 

(If. 19)   ^«A+i  =     ^^""'^^i^yy^^n+i^'^n+i^^n-fi 


C  h2>'  y"(x^^.))  +  (>(h2(l-^)). 


This  is  the  difference  equation  that  would  result  if 
the  initial  value  problem,  the  variational  equation. 


7^ 


(4.20)     n  e'(x)  =  f^(y,x)  e{x)  -  C  [x^  y"(x), 
e(0)   =  0 


were  solved  by  the  finite  difference  methods  under  discus- 
sion. 

The  intial  conditions  for  "e.   are  e".  =  Olh  ~^), 
1   =   0,1, ...,k-l,   where  we  have  made  use  of  (d').   These 
values  differ  from  the  corresponding  values  of  e(x. ) 
(i  =  0,1, ...,k-l)   by  at  most  0{h   ~^).   Hence  by  the 
proof  of  Theorem  4.1  (or  (4.l6))  we  can  write 


e   =  e(x  )  +  (^(h-'-"'*')         for  0  <  n  <_  N 


(4.21)  e^  =  e(x^)h^"^  +  <3"(h^^-^^^)   for  0  <  n  <  N. 


In  order  to  obtain  an  estimate  for  e(x)   we  must  first 
obtain  an  estimate  for  [i     y"(x).   In  Vasil'eva  [4]  it  was 
shown  that  under  the  conditions  of  Theorem  1.2  and  condition 
(c')   y(x,M-)   has  an  asymptotic  expansion  in  \x     of  the  form 


y(x,|i)   ~  <t>(x)  +  9{x,[i)    +    U{[i) 

where      <t>{x)      is   a  smooth  function  and      |0(x,^.)l    <   Ce"  -^^^ . 
Therefore, 
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M-y 


'  =  f(y,x)  =  f(4.(x)  +  9{x,\i),x)   +  ^(n) 


=  f  e{x,[i)  +   (^(ix) 


or 


iM-y'  I   <  C|x  +  C  e 


-Lx/m- 


Similarly, 


(4.22)  ln^y"!   -   |f  uy'  +  M-f  I 

<   Cn  +  C  e-^^/^' 


Integrating  the  variational  equation,  (4.20),  yields 

^^-L{x-t)/kL 
0 


|e{x)|   1  j  ^ |Cn^y"(t)ldt 


<  J  ^ p^ (CM.  +  C  e-^^/^)dt 

0 


Thus 


(4.25)      |e(x)|   <  c^.(l  -  e-L^/^^)  +^e-L^> 


76 


Substituting  (^.25)  into  (4.2l)  gives 


V4.2k)       \e    1   <   C[M.(1  -  e-^^/^)  +^  e'^^/^^]  h^'^'  +  0{h^^^-''h 


Now  (^.24)  will  be  an  improvement  over  {h .l6)    if  l 
is  small,  i.e.,  if  ^  =  1,2,3,  or  ^. 

Remark  k:      If  in  Theorem  4.1  R   eind  T   are 
n        n 

(Xh^P+l^^l-^h   and  e.  =  (:>-(hP^l->'^)   for  i  =  0, 1,  .  . .  ,k-l, 
then  the  same  proof  will  lead  to  the  result  that  e   = 
O'ih^^^-'^h      for  n  =  0,l,...,N. 

We  are  interested  in  investigating  the  "best  method" 
when  the  only  essential  root  is  C,   =  1.      By  the  "best 
method"  we  will  mean  the  one  which  allows  both  ^   and  T  , 
the  truncation  error,  to  be  a  minimum  simultaneously. 

If  the  schemes  with  the  smallest  truncation  errors 
also  have  the  property  that  ^a   ^  0     i  =  0,1, ...k,   then 
by  Corollary  1  Z     will  have  the  value  1.   This  is  pre- 
cisely the  case  when  k  ==  2.   By  using  the  methods  outlined 
in  Chapter  5  on  optimal  methods  we  find  that  the  "best 
methods"  for  the  roots   C  =  1  and  C  =  r  where   ir|  <  1 
take  the  form 


Y  _^o  -  (l+r)Y  ^,  +  r  Y 
n+2   ^     n+1      n 


4?  t{5+r)F^^2  +  (8-8r)F^^^  +  (-5r-l)Fj  +  T^ 
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where 


_   (1+r)   (IV)  A 


Now  if   -1  <  r  <  -  ^  all  p. 's   will  be  greater  than 
—        D  1 

zero.  Therefore,  we  merely  pick  a  value  of  r  close  to 
-1  in  order  to  make  T  small.  Experimentally,  it  was 
found  by  Hull  and  Newbery  [5]  that  r  =  -.97156  gave  the 


smallest  error  for  a  particualr  ordinary  differential 
equation  with  the  property  that   -L  <_  ^  <_  -L  <  0. 

For  the  case   k  =  3   let  the  roots  be   C  =  Ij  ^-^> 
and  rp  with   |r  j  <  1   and   | r^ |  <  1.   After  a  few  pages 
of  calculations  we  find  that  the  optimal  methods  are  charac- 
terized by 


d(C)   =  P^C^  +  ^2^^  +  ^l'^  +  ^o 


where 


P5  =  ^  f9   +  r-L   +  rp   4-  r^rp] 

Pp  =  -^Tj-  [19   -   13(r^   +  rp)    -  5r-^rp] 

^1  "  -^  [-5    -   15(r-j_  +  rp)    +  I9r^rp] 

Po  =  ^  t    1   +  ^1  +  ^2   +  9r^rp] . 
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Clearly  3^  >  0  for  all   |r  |  <  1,   |r  |  <  1.   For  3 
we  observe  that 


19  -  Sr-^rg  >  Ih 


We  must  require  then  that 


1^    -   13{t^  +  T^)      >  0. 


For  example,  if  r,   and  r^   lie  in  the  shaded  region 
of  Figure  1,  then  Po  i  0. 

In  order  to  guarantee  that  P-,  ^  0  we  must  require 
Re  r,  jl  0,   Re  rp  <_  0  and 

-13(r^  +  r^)   >  5 

or 

I9r^r2   >  5  • 

These  conditions  will  be  satisfied  if   r,   and   r^   lie 
in  the  shaded  region  in  Figure  2. 

For  the  case  p   we  notice  rhat  j.f  both  r.,   and  r^ 
lie  in  the  unit  circle  in  the  left  hand  plnne  then  3„  ^  0. 

It  is  clear  then  that  the  condition  that  guarantees 
p,  _>  0  is  the  most  restrictive.   For  example,  the  Adam's 
method   (r,  =  rp  =  0)   does  not  satisfy  this  condition  as 
can  easily  be  seen: 
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X  =  7/13 


Figure   1 


X   =   -5/26 


y  =75719 


y  =  -/57I9 


Figure   2 
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V3  -  V2  =     4?  ^9F^^^  +  19F^^2  -5Pn+l  +  ^n^ 


We  now  investigate  the  truncation  error  for  the  case 
under  study.   Here 


T   =  C,  y(^)  h5 
n      5 


where  C_   is  given  by 
3 


The  minimum  occurs  at 

— ^  -   -19r„  -  11  =  0   or   r,   =   -  ii 


dr^      "^2  ^    ^.    ^2       19 

^  =  -19r^  -  11  =  0   or   r^  =   -  Al 


and 


S 


_11  .  _  J- 


Fortunately,   r,  =  r^  =  -II/19  lies  in  the  region  which 
guarantees  that   P^^O  i=  0,1,2,3.   If  this  were  not  the 
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case,  we  would  have  to  accept  a  larger  value  of  \C^ 


\     ~     ir-^i 


For  the  case   k  =  4   we  know  from  Dahlquist  [ 1]  and 
Henrici  [2]  that   T   is  a  minimum  for  the  roots   C  =  l,- 
C  =  e'^  ,   and  C  =  e""""  .   This  is  clear  since 

^  [_27  -  STr^r^r^  -  ll(r^  +  r^  +  r^ ) 

-  llCr-^rg  +  r^r^  +  r^r^)]  y^^^^x)  h 


where  the  roots  are   C  =  1,  ^-^,    ^2'    ^^^   ^3'   ■'"^  therefore 
seems  obvious  that  we  should  investigate  the  sign  of  P^ 
for  the  case  r.,  =  -s,   rp  =  se   ,  and  r.,  =   se~    where 
0  <  s  <  1.   The  corresponding  P's   for  the  optimal  methods 
are  given  by 


^4  =7^0^251  +  I9(r-j^+r2+r^)  +  llCr^rg+r^r^+r^r^ )  +  I9r-^r2r^] 
P^  =7j^[&46  -  y46{r^+r^+T^)    -  7k{T^T^+r-^r^+T^r^)    -   loer^r^r^] 
^2  "  7io^-264-^56(r^+r2+r^)  +  i[56(r^r2+r^r^+r2r^)  +  264r-^r2r^] 
^1  =  7^0^^°^  "^  T^lr^+r^+r^)  +  ^^eCr^rg+r-^r-^+r^r^ )  -  &46t^t^t^] 
Pq  =7^0^-19  -  iKr^+rg+r^)  -  I9(r^r2+r^r^+r2r^)  -  251r^r2r^]  . 

Now  Pi,  ^  0  for  all  r-j^,  v^,      and  r^  with   jr^l  <  1.   On 
the  other  hand  P^  >  0  if  -1  <  Re  r^  ^  0.   This  is  true 
because 
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>_     &46  -  7'4{j>)    -   106 
>   0. 


For  the  analysis  of  Pp  we  set  r.  =  -s,   r^  =  se   ,   and 


r-,  =  se~    and  ^r  <   9  <  ir 
5  d  —       — 


-26h    +  h^6s{l   -   2COS0)    +  4563^(1   -   2COS0)    -   264s^ 


>   -528   +  ^563   +  ^563^    . 


Let 


h(s)      =     -528  +  k^6s  +  ^56s^ 
h'(s)      =       456  +  912s    . 

Thus  for   0  ;!  s  <_  1,   h(s)   is  rnonotonically  increasing. 
Also  h(.695)  =  0.   Therefore,  the  allowable  values  of  s 
are    .695  <_  s  <  1.   Figure  3  shows  the  acceptable  values 
for  r.  =  -s^   rp  =  se   ,   and  r^  =  se 

p,  >  0  for   -1  <  Re  r.  jf_  0  and  all  r.   chosen  along 
the  same  radii  because 
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Figure  3 


Figure   h 


Q'4 


106  -  7^s(l  -  2cose)  +  3^6s^(l  -  2COS0)  +  6^6s^ 


>  106  -  222s  +  y46s^     =     h(s) 


h'(s)   =   -222  +  692s   =  0 


h(s)   is  monotonically  increasing  if 


222 
s     1     ^     ^      007 


But  h(.307)  >   0.   For   0  <_  s  <  .307   h(s)   is  clearly 
greater  than  zero.   Therefore,   h(s)  >  0  for  0  <  s  <  1. 


Finally,  for  P   our  analysis  is  similar  to  the  above 


where 


•19   +  lls(l   -   2cose)    -   I9s^{l   -   2cose)    +  251s^ 


>      -19   +  lis    -  57s^   +  251s^ 


-76  +  lis   +  251s^      =     h(s) 


h'(s)      =     11  +  753s^ 


3 
Now 


if     s   >^  /  ■^sT  ^    -^72,      h(s)   _>   0.      This   example   is   illus- 
trated in  Figure  4 . 
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We  therefore  see  that  for  the  case  k  =  ^   the  require- 
ment Pp  ^  0  and  3  >^  0   (Figures  3  and  4  respectively) 
provide  the  limiting  choice  of  r, ,  rp,   and  r^ . 

As   k  increases  the  analysis  becomes  much  more  diffi- 
cult and  even  calculating  general  expressions  for  the   P . ' s 

and   T   in  terms  of  r, ,  r„,  ...  ,  r     is  very  tedious, 
n  X   c         K-j. 

We  therefore  resort  to  a  slightly  different  approach. 

From  our  analysis  of  k  =  2,3 >   and   4   we  suspect 
that  for  the   r. 's   in  the  negative  half  plane  and  near  the 
unit  circle  there  is  some  hope  that  for  k  >  4   all   P^^'s 
will  he  greater  than  zero.   We  make  use  of  the  following 
formulas  derived  by  Hull  and  Newbery  [6]  for  optimal  methods 


T   =  R 
^         (k+l)I 


(k+2)  ,k+2     ^  ,  ,^ 


k       i 

R   =  ^   a.    -^  x(x-l)  •  •  •  (x-k)dx 

i-1 


k        1 

i=l  ^  ^  ^. 


where  a.  ,  =  a.  +  a.  ,.  +  •  •  •  +  ot   and 
1-1    1    1+1  k 

p.  =     t~a.         r^  x(x-l)--.(x-k)  ^^ 

for  j  =  0,1, . . .,k  . 


The  above  integrals  can  be  calculated  exactly  by  using 
the  Newton-Cotes  formulas.   We  have  programmed  the  CDC6600 
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computer  to  calculate  the  values  of  p .   for  the  limiting 

J 

values  of  r.   where  we  suspect  favorable  results  for  the 

p.'s;   that  is,  for  k  even  let  one  root  he  at   C  =  1, 
J 

another  at   C  =  -Ij   and  all  the  remaining  ones  at  C,   =  ±±. 

For  k  odd  let  one  root  be  at   C  =  1   and  all  the 
others  at   C  =  ±i-   We  refer  to  this  choice  of  the  roots 
as  a-min.   This  is  in  contrast  to  a-max  where  one  root 
is  chosen  at   C  =  +1   and  the  remaining  ones  at   C  =  -1- 

Although  both  a-min  and  a-max  define  unstable  schemes, 
in  practice  we  would  choose  root  at   C  =  1   and  the  other 
roots  inside  the  unit  circle  but  near  the  roots  of  a-min 
or  a-max  when  they  lead  to  p.  's  >^  0. 

The  results  of  the  computer  calculations  are  given 

in  Table  1  for  k  =  3,  ^  ,  .  .  . ,  l8.   By  P  _>  0  we  mean  that 

^ .    >_  0     for   i  -   0,1,..., k.   P  <  0  means  "Chat  at  least 

one  p.   was  less  than  zero. 
1 

We  will  conclude  this  section  by  saying  that  at  least 

for  the  case  k  even  we  would  expect  T   to  be  small  for 

10^     ±9  i0(k-2)/2 

L,   =   -s,   se   ■",  se    ,  .  .  . ,  se   ^    '    where   s  <  1,  but 

near  1.   For  the  cases  k  =^  2,3,  •••j7   if  we  choose 


27rim./k         ^  ^     ,  ., 
^m     ^     ^^  ^  ^  l,2,...,k-l 


the   p.  's  _>  0  for   s   sufficiently  large  because  for   s  =  1 
these  roots  generate  the  Newton-Cotes  formulas  where  the 
p.  's  are  positive  for  k  <  8. 
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Degree  ct   -  max  a   -  min 


k 

= 

3 

P 

> 

0 

p 

> 

0 

k 

= 

h 

P 

> 

0 

p 

> 

0 

k 

= 

5 

P 

> 

0 

p 

> 

0 

k 

= 

6 

P 

> 

0 

p 

> 

0 

k 

= 

7 

3 

> 

0 

p 

> 

0 

k 

= 

8 

P 

> 

0 

p 

> 

0 

k 

- 

9 

6 

> 

0 

p 

> 

0 

k 

= 

10 

P 

> 

0 

p 

< 

0 

k 

= 

11 

P 

> 

0 

p 

> 

0 

k 

= 

12 

P 

> 

0 

p 

< 

0 

k 

= 

15 

P 

> 

0 

p 

< 

0 

k 

- 

14 

P 

> 

0 

p 

< 

0 

k 

- 

15 

P 

> 

0 

p 

< 

0 

k 

= 

16 

P 

> 

0 

p 

< 

0 

k 

= 

17 

P 

> 

0 

p 

< 

0 

k 

= 

18 

P 

> 

0 

p 

< 

0 

Table 

1 
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CHAPTER  V.   EXTENSION  TO  HIGHER  DIMENSIONS 


Consider  i-L  -^  =  f(y,x)   where 


y(x)  = 


y^(x) 


y  (x) 


and  f(y,x) 


a,   1    2         S    ^ 

f  (y  .y  ,  ■  ■  ' ,y   ,x) 

^2/12       s   X 

f  (y  ,Y  ,  •  • • ,y   .x) 


„s,  12      s   N 

f  (y  >y  ,••-,¥   ,x) 


¥e  cam  use  the  same  notation  as  previously  provided  we 
make  the  following  definitions: 


Sf^(y,x)   af^(y,x) 


fy(y,x) 


^y' 


^y 


2 


aff(y,x)^  Bf^(y,x) 

Sy      Sy 


af^(y,x)   Sf^(y,x) 


^y' 


^r 


Mi(y,x) 
^y' 

£i_(y,x) 


af°(y,x) 
3y" 


Similar  meanings  will  be  prescribed  to  f   ,   f   ,   and 

XX    xy 

f   .   Also,  we  have 
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n 


y^(xJ 


n 


y^(x„) 


y^n) 


n 


n 
n 


n 


n 


f-'-(Y-^,Y^,  ...,Y^,x    ) 


f^(Y-^,Y^, 


,Y^,x    ) 
'    n'    n 


n 


f(y(x^),x^). 


<l>(x^)      = 


.l(xj 


.2(xJ 


*    (x    ) 
^    n 


and        f 


■yn 


■y|x=x 


n 


Finally,  we  replace  the  absolute  value  with  the  maximum  norm: 


n 


max    I  yj;  | 
i=l,2,...,s  " 


y 


n 


Let     A  =   (a.  .)      be   aji      s     by     s     matrix.      Then 


A 


A 


=     max 


a 


i     k=l 


ik' 
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Using  this  notation  we  consider  the  linear  system  defined 
by 


^  =  f(y,x)  = 


1  ^     2  ^ 


a2^y  +  a^gy  + 


dx 


■^  ^is^^ 


+  a_  v"" 
2s"^ 


^sl^   +  ^32^   + 


+  a  y 
ss*' 


or  ii-T^  =  Ay  where  A  =  (a.  .)   is  a  constant  matrix 
dx    "^  ij 

The  related  difference  equation  is 


(5.1)    a,  Y  ,,     +  a,  ,  Y  , .  .  +  •  •  •  +  a  Y 
^•^       k  n+k    k-1  n+K-1  o  n 


1  _ 


=  h-^(aAY^^,  +  Pv_,AY^^._,  +  ...  +  p  AY  ) 


k  n+k   ""^k-l  n+k-1 


o  n 


Let   P  be  the  matrix  which  reduces  A   to  its  Jordan 

Canonicax  form.   Let;   PV   =  Y   and   P~  AP  =  J.   It  should 

n    n 

be  noted  that   P   is  a  constant  matrix.   Substituting  into 
equation  (5.I)  and  multiplying  on  the  left  by  P~  ,   we  get 


k  n+k    k-1  n+K-1  o  n 


=  hi->'o^jv^^^  +  \_iJV^^i,_i  +  •••  +  Po^J 
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First  of  all  let  us  assiome   J  is  diagonal  with  eigenvalues 
A^   (i  =  1,2,.  ..,s). 

We  look  at  the  component  form  of  the  last  expression: 


a.  V^^,  +  .  .  .  +  a  V^   =  h  ^^O,  A.V^^,  +  •  •  •  +  p  A.V^) 
k  n+k  on         ^  k  i  n+k         ^o  i  n 


By  Theorem  2.1  for  h   sufficiently  small  the  roots  of 
p(C)  -  h-^"^A.6(C)  =  0  are  of  the  form 

-.         1/m. 
(5.2)  Cj   =   Cj^  +  ^jPr~ +  (7(h2^1-^^/^0) 

where   p(C-  )  =  0  and  m.   is  the  multiplicity  of   s-  • 

^^  JO  J  ^      "^        JO 

However,  we  have  imposed  no  conditions  on  A.   as  yet 

Since   C-   =1   is  always  one  of  the  roots  of  p(C)^ 

we  must  first  stipulate  that  Re  a .  <  0 .   As  in  the  one 

dimensional  case  if   |C-  I  *=  1,   then 

'  Jo' 


lim  K  .r  =  0 

h  -*  0   ^ 

X  =nh 
n 


Assume  now  C  •   =  e    and  hence  m  =  1.   Let   A. 

JO  X 

-A.   +  iA.,   where   A.   >  0.   Recall  that 

lO       ll  lO 

(5-3)  ^  =  Pi  +  ^^1  • 
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Equstion  (5«2)  may  now  be  written  as 


(5.^)    ^j   =  COS0  +  (A^^Pj  +  Aiiqj)h^"^ 


+  i[sin0  +  (A.^q   -  A.^p  )h^->']  +  Olh^^^-'^h 


and  consequently 

(5.5)      Kj.|^   =   1  +  [2(A^^Pj  +  A^^q^.)cos0 

+  2(A^^q^.  -  Aj_^p^)sin0]h^~^  +  ...  . 

In  order  that   lim  |C-i    converge,  we  must  require  that; 

h  -*■  0   ^ 

X  =nh 
n 

the  term  in  brackets  in  (5.5)  be  negative.   In  general  this 
is  not  a  very  practical  condition,  since  it  depends  upon 
the  eigenvalues  of  A,   which  are  difficult  quantities  to 
obtain. 

We  may  obtain  a  more  practical  condition  by  requiring 
that  the  eigenvalues  be  real  and  negative   ( A . .,  =  O)   when 
we  wish  to  consider  difference  schemes  having  roots  of 
p(0   on  the  unit  circle,  but  not  identically  equal  to 
unity. 

When  A   =0   (i=l,2,...,s)   we  have 


(5.6)  ICj-l^   =  1  +  [2(A^^PjCos0  +  A^^q^.sin0)]h^"^  + 
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Our  condition  now  'becomes 


(5.7)  p  .COS0  +  q.sine   <   0; 

J         J 


that  is,  we  require  (X-stabllity  in  order  to  guarantee 
convergence . 

The  following  theorem  will  be  useful  in  order  to  assure 
that  the  eigenvalues  of  A  are  real  and  negative: 

THEOREM:   Suppose  that  A  is  syrnmetric,   a^^  <  0  and 

s  , 

a..  +  2—    la.  .  I   <   -M   <   0    i  =  l,2,...,s 
11    j=l   ^^ 

where  the  prime  on  the  summation  sign  denotes   i^iil   ^^ 
not  included  in  the  sum.   Then  the  eigenvalues  of  A   are 
real  and  strictly  negative. 

Proof:   See,  for  example,  Varga:  Matrix  Iterative 
Analysis . 

We  give  an  example  of  a  scheme  which  is  ^.-stable  but 
which  does  not  give  convergent  results  because  the  term  in 
brackets  in  (5.5)  for  this  example  is  positive,  i.e.,  this 
method  can  be  used  for  the  scalar  boundary  layer  equation, 
but  not  the  vector  one . 

Example  1.   Let   p(  O  =  C"^  -  C^  +  C  -  1-   The  roots  are 
C  =  1,  ±i,   and  the  M'-stability  condition  implies  that 


9'-i 


1  >  Pi  4-  Pg^ 


Let   d(C)  =  I  C^  +  I  C  +  1.   Therefore, 


-6(1) 
P'CD 


1   1  . 
IT  -  IT  ^ 


Suppose  the  Jordan  Canonical  form  of  A  Is 


-2-'4±  0 

0     -2-f4i 


Substituting  the  eigenvalue   -2-^i  and   C 


JO 


i  into  (5-5)> 


we  obtain 


In  order  to  prove  the  analogy  of  Theorem  ^ . 1  for  the 
vector  equation  M-y'  =  f(y,x)   we  have  seen  that  in  certain 
cases  it  is  necessary  to  restrict  the  class  of  acceptable 
functions  f(y,x) .   First,  as  before,  we  assume  f  (y(x),x) 
is  continuous  and  bounded  in  D.   We  distinguish  between 
two  cases: 


(5.8)   Case  1.   The  roots  of  p(C)   lie  inside  the  unit 
circle  with  the  exception  of  C  =  1-   Here  we  require  that 
the  real  parts  of  the  roots  of  the  characteristic  equation 
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Det  11  ^(y(x,^^),x)  _ 
"  oy 


=  0 


are  negative  and  bounded  away  from  zero  in  D  for  all 
l-L   sufficiently  small. 

(5-9)   Case  2.   The  roots  of  p(C)   other  than  C  =  1 
may  lie  on  the  \anit  circle.   Here  in  addition  to  [X-stability 
we  must  require  that  the  roots  of  the  characteristic  equa- 
tion are  strictly  negative  and  hoimded  away  from  zero  for 
all  values  of  x   in  D  and  |J.   sufficiently  small. 


Let  us  now  investigate  the  form  of  Duhamel's  principle 
for  the  consistent  difference  equation  satisfying  condition 
(5.8)  or  (5-9) 


k 


k 


(5.10) 


.^-^  E  Pi' 


a.Y  ^-  -  h   ''  7   P.JY  ^.   =  Q  , 


J  is  a  matrix  in  Jordan  Canonical  form.   We  claim  that  the 
solution  to  (5.10)  takes  the  form 

r   n-k 


(5.11) 


n 


/   %      .  ,  Q.  +  y      n  >  k 


J-0 


■J-1   J    n 


1 


n 


n  <  k 


where  ^   is  defined  as  a  matrix  satisfying  the  following 
equation: 
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0 


(5.12)  Ylln+±^^±   -   h^'^V^   = 


i-0 


with  initial  conditions   ©^  =  ^-,  =  •  •  •  ^,  „  =  0   and 

-*-o   -"^l       -'^k-2 

^,  _-,   is  chosen  to  satisfy  the  equation 


1-7, 


(5-15)   ^i,_i(a^  -  h-^-V)   -   !        ■  •  .     i   =   1 

°  1 


Since   J   is  in  canonical  form  we  see  immediately  that 

Ik-1   =   0   if   i  >  j 

(5.1^) 

i-i  =  («k  -  ^''Vk)-'- 

The  other  elements  of  the  matrix  ^,  -,   will  depend 
upon  the  particular  form  of  J. 

f        is  the  solution  to  the  homogeneous  equation 
((5. 10)  with  0^  =  0)   with  initial  data  Y  ^Y^-'-^-^.i- 
Therefore, 

k-1  k-^-1 
(5.15)  \     =  5  '  5  '""^-J  "  '^'"V/'5n+fc.i.j.i)Y, 

for  n  >  k 


where  ^   is  defined  by 
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k_ 


(5.16)  Yl  («i  -  h^"^PiJ)^+i 


with  initial  conditions  ^  =  !]_  =  •  •  •  =  $^_2  "  *^  ^^^ 
¥,  -,   is  chosen  to  satisfy  the  equation 

(5.17)  («!,  -  h^"X"^^^k-l   "   ^• 

To  verify  that  (5-ll)  is  really  the  solution  to  (5-10) 
we  consider  the  homogeneous  solution,   ^  ,   first: 


7    (a.  -  h^-^p  J)^ 
i=0    ^ 

k  k-1  k-^-1 


i=o   "-       ^   ^=0   j=0   "^  J        "^  ^ 


5i+i+k-^-j-l^^^ 
k-1  k-^-1  k 

H  (  H  ZI  («k-j  -  ^^"X-/^("i  -  h'"''PiJ) 

^n+i+k-^-j-l^^i 


0     for     n  f^  0 
k 


and   0<^+j<k-l. 


hecaus    

i=0 
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For  n  <  k,   we  have 
k-1  k->g-l 

} 


IS.-X    J\.-£-X 


1-7 


^  k        k  -^k-l  n 


=  \ 


On  the  other  hand,  the  particular  solution  satisfies 
the  equation  (5.IO)  because 

n-k  n-k        k 

ZZln-j-A  =       Uln-j-l^  ^  (^'i  -  h^'^'^iJ^Y  i) 


Let   j  +  i  =  ^ .   Then 
n-k  k  n+i-k 


V-7T  -*-n-.-]-l  .1     ^: — K-  V-r  -"-n-^-l+i   1        1 


j^     --j-1^0     ^^^n-^-l+i^  i       ^i"'^^ 


But  for  the  particular  solution  we  impose  the  initial 

values   Y^  =  0   for   i  <  k.   Since  initially  ^  =  ^-[  =  •  •  ' 

=  $,  „  =  0,   we  have 
-'-k-2 


n-i^.  K   n-M 
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k    n 
i=0   i=0 


Zi  ZI^n-^+i-l^°^i  -  h'-^PiJ)Y, 


where  the  last  result  follows  from  the  definition  of 
^    =0  for  n  =  1,  2,  ...  . 
Hence, 


n-k 


n    k 


j=0 


^n-j-l^j   - 


JZ    IIln-^+i_i(«i  -hl"V^Y^ 


k 


1-7. 


i=0 


^  ,  ^.  (a.  -  h  ^^p.  J)Y 


1-7. 


=   ^K-l^^k  -  ^   'V)\ 


n 


hy  ( 5 • 12 )  and  { 5 • 13 ) • 


It  is  now  necessary  to  prove  lemmas  analogous  to 
Lemmas  1  and  2 .   We  will  merely  illustrate  the  proof  for  a 
special  case;  namely. 


(5.18) 


J  = 


\ 

1 

0 

0 

0 

\ 

1 

0 

0 

0 

\ 

0 

0 

0 

0 

"-2 

The  more  general  Jordan  form  lends  itself  to  a  similar  proof. 
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LEMMA   1'.      Let   the   consistent   difference   equation 
k 


(5.19)  ZIln+i^°^i   -   h^'X^^ 


0 
i=0 


satisfy  condition  (5.8)  or  (5-9)  and  initial  conditions 

lo  =  ll  =  •  •  •  =  lk-2  =  0  ^^  lk-1  =   (^'k  -  h'" V^"""- 
Then 


N 


7~  llf  I   1  ^     J,k  =  1,  2,  ...  ,  s 


1^   "       h-^'^ 


for  h  sufficiently  small  and  N  any  integer. 

Proof:  We  supply  the  proof  for  the  case  when  J  has 
the  form  given  by  (5.18).  The  initial  conditions  on  ^■^_-^ 
are  given  by 


$k-i'°k  -  '''''V   =  1  • 


Equating  terms  gives 


i^-i  -  ^\  -  ^'-Wi^-'     ¥A  -  W-i '^  t^-i  - 


0 


^k- 


^^-''"'■^  i.fx  =  (".-^-Vi'"  i.'-M.:?=° 


1   (a^-hl-'-A.P^)^ 
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lc!i  =  (v^'-Vi'"'     i.-'i  =  ° 


.^l-i,?l  =  1^1  =  0     i.*i  =  K  - -^^-v.'-^ 


With  these  initial  conditions  it  is  clear  that 


J  ^J-  =  0   for   i  >  J  and  n  =  0,  1,  2, 


Also,   $  i^  =  I  ^^  =  I  ^^  =  0   n  =  0,  1,  2,  ...  . 


By  Lemma  1  and  conditions  {5-8)  or  {ti-9)    on   A-   and 
Ap   it  follows  that 


n-1 

1^  JJl   ^    C 


< 


i^o  "^  ^   -  h^-^ 


j  -  1,2, ... ,s 


The  remaining  elements  of  ^  satisfy  the  following 
equations : 


a) 

k 

>        («i 

i=0         ^ 

-•^'"Vi'inii  = 

=  hi-'-)    p.fH 

b) 

k 

-   hl-rp^A,)|^g      = 

=      hi-'')        Pjl2 
i=0          n+i 

c) 

k 

-h^-Vi)|.?i  ^ 

=  >^"^>1paS 
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We  solve   a)  by  Duhamel's  principle 
n-k  k 


-^   n      ^-TTT^  -J-n-j-l      4—^  i-*^,]+i     n      — 
(5.20) 


j=0  ^^^J^-^      i=T) 


$  n  <  k 

n 


where  ^   the  one  dimensional  Green's  function  defined 
in  Lemma  1.   We  need  the  following  theorem  from  advanced 
calculus: 

CO 

THEOREM.   Assume  that   / ^   converges  absolutely 

n=0  ^ 

and  has  sum  A,   and  suppose         / b   converges 

00  n=0    CO 

with  sum  B.   Then  /   c   where   c^   =   }    a,b.^  - 

feo  ^        ^     fco  ^  ^-^ 
converges  and  has  sum  AB. 

Proof:   See  Mathematical  Analysis  by  Apostol. 
From  (5-20)  we  have 

n-k 

A 


(5.21) 


j=o   "■"^' 


n-1 


J=0 


IX  11 


since   |^  .  |   Is  decaying 


N  N  N 

fc^  '^-'  -     h^->'  h^-^   feo  '  -'  -  h^-^  U  -' 
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A 

where  we  have  made  use  of  Lemma  1  on  the  functions  ^^ 

and   Q    and  the  above  stated  theorem.   Now 
^n 


A 
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n     ^  n 


=  c^^'.r     n  <  k 


and 


i  =   !l  ('|l'(°K-j  -  •^'-^p.-jH'U-i-j-i'if 


,  A 

=  ^ ^  n  ^  k 

— 1  p 
where  we  have  used  the  initial  conditions  for  ^^ 

Therefore, 
N  1  y  N 


<   C      by  Lemma  1  again. 

We  therefore  have  that 

N 

7F12 1   ^    C 


II  lln 


n^  ""   "  h^''' 


10^ 


Equations  b)  and  c)  are  treated  in  the  same  way.   Lemma  1' 
is  now  complete. 

LEMMA  2 ' . 


—  •1  C  -,  -I     n    -, 

llf  I   <  [ r^n^-1  4-0(1  -  h^-^L)  n^-1] 

^    -  ^(l-7)(m-l)/m 

for  all  n  and   j,k=l,2,...,s  . 

Proof:   For   l^""""^!   the  proof  is  the  same  as  Lemma  2. 
From  (5-21)  and  Lemma  2 

n-k  ^ 

Ifl^l   <  h^-^p  y~  1$   .  J(lF.^l+---  +  lK^|)  +  1^  I 
iXn  '   —      ^   4^  '-^n-j-1'  ^  '-^  J  '      '-^j+k'     '  n' 

n-k 

<   (k+l)h^->'p  ZI  t ^^""^   +  C(l-h^-^L)^-J] 

j-0  ,  l-7(m-l)/m 


( Cr^ +  Cd-h^-^L)"^)  +  1$^ 

^(l-7)(m-l)/m  ^ 


n  ,    n 


<  2_^^^ +nC(l-hl-^L) 


^(l-7)(m-l)/m 


where  we  recall  that  ^  -     =1 — -^^ —  ¥ 

t15 


For  ^  -^  we  have 
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\t^\     1     hl-^P^I  lin.j.il(llfl---lif.j)    -   l^n! 


2_,   n  „  -,  n 


^l-7(in-l)/in 


LEMMA   2" 


^^  ^    -     ^(l-7)(m-l)/m 


for   all     n      j,k-l,2,...,s,      0<r<l,      L>0      and 
m  =  max  ( m-,  ,  irip ,  .  .  . ,  m . )  . 

Proof:   The  proof  follows  the  same  kind  of  arg-oment 
as  the  ones  found  in  Lemma  1'  and  Lemma  2'. 

With  the  aid  of  these  lemraas  we  can  prove  the  following 
theorem: 

THEOREM  5-1.   Let  the  consistent  finite  difference 
equation 

k  k 

J        a.  Y  ,.    =    h^"^  /    p.  A  Y^  ,,  +  R 
A — 7:r  1   n+x  4 — n   1    n+i    n 

1 =0  1 =0 


where  A  =  (a.  .)   is  a  constant  matrix,  satisfy  the  following 

-^  J 

conditions: 


io6 


a)  R^   -  (5^(h^'^^)   as  h  -^0 

b)  Condition  (5>8)  or  condition  (5-9)  is  satisfied 

c)  |Y.  -  y(x.)|   <_  Th-^"'^  for  h  sufficiently  small 
and  i  =  0,l,...,k-l.   T  is  a  positive  constant  indepen- 
dent of  h. 

Then  the  difference  equation  is  convergent. 

Proof:   Subtracting  the  difference  equation  satisfied 

by  y   from  the  above  and  setting  e   =  Y  -  y   gives 
•^   •'n  '^   n    n   "^n  '^ 

k  k 

y^  a.e  ^.   =  h^"^  y~~  p.  A  e  ^.  +  R   -  T  . 
t^  1  n+i  f^   ^1    n+i    n    n 

If  we  let  e   =  Pu   where  P~  AP  =  J  and  multipliy  on 
the  left  by  P~  ,   we  obtain 


i~0   -^  ^^'^  i=o 


^   a.u  ^.   =  h-^"^  )   P.  J  u  ^.  +  P"^(R„  -  T  ) 
4 — TT     1  n+i  4 — TT     1    n+i      ^  n    n 


The  proof  now  follows  the  same  course  as  Theorem  2.2  with 

the  use  of  Lemmas  1',  2'  and  2"  instead  of  Lemmas  1  and  2 

/^n  h             /^n  h 
However,  we  must  choose  N  =  s  n  instead  of  n       

because  of  the  slight  difference  between  Lemma  2"  and 

Lemma  2.   This  becomes  clear  when  we  note  that 


,  s-1  n 

^s-1  n        x_  r       ^  -I 

_n I <  J2 <  x^"^  <   1 

(l-7)(m-l)/m  ~    h^  ~   "^    ~ 


h 


and 
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lim  h^  ^  ^n  h   =   lim  -^  h^~'^      =     0. 


h  -^  0 


h  —  0 


7-1 


Q.E.D. 

THEOREM  5.2.   Let  the  consistent  finite  difference 
equation 


k  k 

(5.22)    II  a  Y     =  h^-^H^ 
i=0  -^  "  -^  i=0   ^ 


F  ^  .  +  R 
n+i    n 


where 


F 


n 


f'(Y^,Y2,...,Y^,x^)  ;| 


f^(Yl,Y2, 


,Y^,x  ) 
n  n 


satisfy  the  following  conditions: 

a)   R^  =  C?(h^^-^"^^)   as  h  -*  0 

Sf   hf 


n 


^2 
b)   The  matrices  ^,   ^,   and  — p  are  continuous  and 


S/ 


bounded  for  0  Jl  x  <_  1  and  -oo<y-L<oo   i=l,2,...,s 

c)  Condition  (5«8)  or  condition  (5 -9)  is  satisfied. 

d)  le.  I  _;^  Th  "■^   for  h  sufficiently  small  and  i  = 
0,1, ...,k-l.   T   is  a  positive  constant. 

Then  there  exists  a  constant   C   such  that   |e  |  _^  Ch  "-^   for 


n  =  0,1,...,  N  where  x^.  =  1> 
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Proof:   The  exact  solution  of  [xy '  =  f(y,x)  satis- 
fies the  difference  equation 

k  k 

(5-23)    /       ct.y  _^.  =     h^"^  }   p.f  ,.  +T 
^         4 — pr  i"^n+i  4 — pr  1  n+1    n 

1=0  1=0 

where  T   can  be  shown  to  be  Cy(h  ~  '^)   as  h —*  0  in  the 
sajne  way  as  in  Theorem  '4.1. 

Subtracting  (5-23)  and  (5.22)  and  letting   e.  = 


Y.  -  y. ,      we  obtain 
1   "^i' 


k 


a.e  ^.   =   h-^"^  )    P.(F^.  -f^.)+R   -T 
.  _^      1  n+i  4-^  "^i'  n+i    n+i     n    n 


In  component  form  we  can  write 

k  k 

r^a.  e  ^   =  h^-^^P  (Ff  -  f  f  )  +R^  -  T^ 
v^  1  n+i         4-rQ  i  n+i    n+i     n    n 


for  £   =   1,2, . . . , s 


k      s  -s-^X 
i^-^ZIp-  II -r'  e  J.  4-R^  -  T^ 


where 


— ^ 

^f   ..  >^-2,   1    ,q1-2    1           S    ,aS£        S                         \ 

n+1  of  (y  ,.+0  ,.e  ,.,..., y  ,.+0   .e  ,.,x  ,.) 

— ^.  =  — H-^^n+i  n+1  n+i'    '"^n+i  n+i  n+i'  n+i 


SL      =     1,2, ...,s 
J   ~   1,2, ...,S 
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Therefore,  in  vector  notation  we  have 


(5-24) 


i=0 


a.  e  ,  . 
1  n+i 


=  h 


'-'  E  Pi 


i=0 


hf      . e  , . 
-s— n+i  n+i 


R   -  T 
n    n 


With  the  choice  of  N  =  s  -J^ —  we  can  find  a  bound  for 

o     ^n  r 

le  I   for  0  <  n  <  N   as  was  done  in  Theorem  4.1.   Assuming 
'  n         —   —  o 

this  has  "been  accomplished,  we  add 


k 


o  n+i 


i=0 


i=0 


to  both  sides  of  (5.24)  and  make  the  substitution  e   = 

Pu   where   P~  L  P  =  J,   the  Jordan  Canonical  form  of  L  ■ 
no  o 

After  multiplication  by  P~  ,   we  have 


i^ 


(5-25)  ZZ  («i   -  h^'^P^J) 


u 


n+i 


h 


1-7 


i=0 


a     •n~l/  Sf      ,   • 
^i^         (^+^ 


5^o)Pu   ^.    +  P"^(R 
oy  n+i  ^    n 


Tn> 


We  let  0   denote  the  right  side  of  (5-25) • 

Using  condition  b)  of  the  hypothesis  we  can  write 


^f^ 

?)f^ 

n+1 

o 

ayJ 

SyJ 

af^^.     Sf^ 

I n+j.       o 

ayJ        ay^' 


-(?(K,,\) 


<_     s+a(|u„,J) 
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where  the   e  can  be  chosen  as  small  as  necessary  by  con- 
tinuity. 

A  boiond  for  Q   is  then 

k 


iQn'   1  h^"^  H  ip-^llPl(lP.lse) 


i=0 


n+i 


k 


+  ^(h^^l-^))  +^(hl-^ZI 


i=0 


n+1 ' 


where  0  {h   ^    "'^  )   is  due  to  the  round-off  and  truncation 
errors . 

The  solution  to  the  system  (5-25)  was  given  earlier  as 


r  n-k 

)   f   .  .  Q.  +     "^ 
h-rr  ^n-,1-1  ^,1      n 


(5.26)     u 


j=0 


n 


^ 


n 


n  >  k 


n  <  k 


If  we  now  introduce  the  max  norm  U  =  msix 

n    .    -, 

1  =0 , 1 ,  .  .  • ,  n 


Uj_  I   and 


follow  the  same  procedure  as  in  Theorem  ^.1  but  using 
Lemmas  1',  2',  and  2"  instead  of  Lemmas  1  and  2,  we  obtain 
the  same  basic  inequality 


K^U^   K^h"^"^   K,,Th^  ^ 


n  —  1-a     1-a 


+ 


1-a 


0  <_  n  <  N 


where  N,   is  determined  by  the  choice  of  e  and  continuity 
considerations . 


Ill 


The  proof  is  completed  exactly  like  Theorem  ^1.1.   Q.E.D 
Since  the  estimates  used  in  Theorem  5.2  are  rather 
complex,  it  is  not  a  simple  matter  to  estimate  the  n^uuTiber 
of  steps  i      occuring  in  Theorem  5 -2.   Instead  we  show  hov; 
better  error  bounds  can  be  obtained  when  we  consider  a  very- 
restricted  class  of  difference  schemes. 

THEOREM  5.3-   Let  the  consistent  explicit  difference 
equation 

k  k-1 

4^  1  n+i  ^  ^1  n+i      n 


satisfy  the  following  conditions: 

a)  IR   -  T  I   <   K  h(P+lHl-7) 
'  n    n'   —   1 

b)  e^   =  ($^(hP^^"^^)   for   i  =  0, 1,  .  .  .  ,k-l 

c)  The  matrix  -r—  is  continuous,  bounded  and  strictly 

negatively  diagonally  dominate  for  all  points  of  interest, 

i.e.,  if  A  =  (a..),   then 

-'-J 


a..   +7    |a..i    <    0     for  all   i. 


H 


1 


d)   a^  =  1,  -a^_^,  -a^_2,...,-a^,  -a^   >   0   and  ^^_.^, 

Pk-2'"-'Po  ^0- 

Then  the  difference  equation  is  convergent. 
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Proof:   If  we  subtract  the  difference  equation  satis- 
fied by  y   from  the  above  and  set  e  =  Y  -  y^-,>   ^e 
obtain 


k-1 


'n+k 


/    [  -a .  e  ,  .  +  h-^"^p  .(F,.-f,.)]+R   -T 
4 — ■pr  1  n+1        1  n+i    n+i      n    n 


i=0 


k-1 


=  II[-a.  +hl->' 


Sf 


n+i- 


i=0 


loy     '^■^^  '^ 


'n+i    n    n 


where 


Sf 


^' 


■n+i  = 


— -,  ^  n+i '  n+i  , 


S,  rS 


af^(C 


n+i '  n+i  , 


— _  ^  n+i  n+i 
ay^ 


^'(cL..x_j  , 


n+i  n+i ' 


oy 


By  condition  c)  there  exists  a  number  M  >  0  such  that 
for  all  numbers   a  >  0  and  b  >_  0  and  h  sufficiently 
small 


al  +  bh^"^  l^+il  < 
dy 


a  -  b  M  1'}'^ 


for   n  =  0,1, ... ,N 
and   i  =  0,1, ... ,k 
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Taking  norms  and  using  the  last  inequality  we  obtain 
k-1 
i=0 


e„,,l  1  II   (-i  -  PiMhl-r)|e„,J  .   K^h(P«)(l-7) 


where  conditions  a)  amd  d)  were  used 
We  now  prove  by  induction  that 


K.hP^l-^) 

1   <  ^ 

'n+k'   -    k-1 

M 


^h 


i=0 


We  assume  this  last  inequality  is  satisfied  for   l®o'' 
|e,  1  ,  .  .  . ,  I  e,  -I  I  .   If  it  is  not,  choose  a  larger  K-,  .   By 
induction  then  we  assume  the  inequality  holds  for  \^^  > 

l^k+l'-'-'^n+k-ll'   *^^^ 

k-1  ^  uP(l-7) 

1=U  V 

i=0   ^ 
„  , p(l-7) 


M  

i=0 


K  h 


-•=0   ^ 

Pd-r) 


i=o  ^ 


11^ 


where  we  have  used  the  consistency  condition  / (-oc.  )  =  1. 

Thus  the  proof  is  complete  by  induction. 


k-1 

H 

i=0 
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VI.   BOUNDARY  LAYER  EQUATIONS  OF  HIGHER  ORDER 


Suppose  we  wish  to  consider  the  second  order  equation 


2 
(6.1)         M.  ^  +  a{x)  II  +  t)(x)  z   =  r{x) 

dx 


Let  g  =  y. 


Equation  (6.1)  is  equivalent  to  the  system 


M-  -^  =  -a(x)  y  -  t)(x)  z  +  r(x) 


dx     ^ 


We  are  therefore  led  in  a  natural  way  to  investigate 
the  system 


M.  g  =  f(y,z,x) 
(6.2)  0  1  X   <   1 

§1  =  s(y,z,x) 


with  initial  conditions 
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y(o) 


y 


o 


z(0)   =   z° 

Following  our  previous  approach,  we  first  consider 
the  linear  case 


M- 


dZ  = 


(6.3) 


dx 


dz 
dx 


a-^-j^y  +  a-^gZ 


ag^y  +  a^gZ 


where   (a.  .)   is  a  constant  matrix.   From  the  ODE  theory 
we  must  require  that  -^  =  a, ,  <      0 . 

Dividing  the  first  of  the  above  two  equations  by  m- 
and  using  vector  notation,  we  obtain 


dy 

dx 

li  y 

dz 

=      A       " 

dx 

II 

where 


V  " 


^11   ^12 


^21  ^22 


a 


11    ^12 


7 


7 


21    22 


s  A, 


h 
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The  characteristic  equation  of  A   is 


a 


11 
p. 


(^  -  >^)(aoo  -  ^)  -  3-21  !l2   ^  Q 


22 


M- 


and 


(6.5)   A  = 


2il 


where 


A   = 


i    a-j_^    a^2 


^21    ^22 


(6.6) 


a 


1      IJ. 


11 


as   li.  -►  0 


(6.7)  A, 


^La^g  +  a-^-^  -[M.a22  +  a^i  -  2M.  det  A(M-a22  +  a^.^) 


-1 


2M. 


^(4[i  det  A)^(M-a22  +  a^^)  ^  +  ...  ] 
_  211 


(6.8) 


det  A  ■      ,,  _  n 

A^  ~  as   M-  -^  0. 

2      a^-^ 


We  note  that  for  M-   sufficiently  small  the  eigenvalues 


are  distinct  and  Ag   is  bounded. 
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The   eigenvectors   are   found   as   follows; 


A. 


1^ 


A. 


'1 


\x      ^1  M-      ^2 


\h 


^21^1      "*"  ^22    ^2      ~      ^1^2 


^21^1   "^  f'22^2      "      ^2^ 


M-agg   +  a^^ 


+  y{[ia^^  +  a^^)^   -  ^ixdet  A_ 


2[i 


^2 

T7 


2m.  a, 


21 


■^  [a^^^   -  txa^g   +y(iia22   +  a^-^^)      -   ^p.   det  A] 


Letting      ^,    =  1     we   obtain  the   eigenvector 


1  ■ 

in 


Remark:   lim  ^„  =  0.   Hencej  the  eigenvector  is 

bounded  for  M-   sufficiently  small. 

The  eigenvector  corresponding  to  A^  is  found  similarly: 


^•11      ^12 


w 


^21  ^1  +  ^22  \      =     W 


^11      ^12 


r]^[\ia,^^  +   a^^  -  /(ixa^g  +  a-n^^-  ^^^det  A] 


2M- 
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-2a 


12 


^n  -  ^^^oo  -^y^^^no    +  a-n)^  "  ^t^  det  A 


^11   ""^22  ■  ^  ^'-"•22  '  ^11' 


li  h 
Letting  t)  =  1  we  get  the  eigenvector   | 

i  1 


,   and  we 


note  that  for  M-   sufficiently  small  this  vector  is  bounded. 
Let 


(6.9) 

then 

(6.10) 


P  = 


1  ^- 
^2   l' 


-1 


1  -  \^2 


1    -11. 
-^2    1 


Clearly,  for  M-   sufficiently  small  P~   is  bounded. 
Also, 


(6.11) 


.-1 


V 


P   =  D   = 


A^    0 
0     A. 


for  M-  sufficiently  small. 

The  difference  equation  associated  with  system  (6.^) 
is 


> 


^n-t-i 


n+i 


k 


=  h 


n^. 


A 


i=0 


h 


n+i 


n+i 


l20 


Let 


Y 


n+i 


n+i 


;!  U 


n+i 


V. 


n+i 


Substituting  this  expression  into  the  difference  equation 
axid  multiplying  on  the  left  by  P   ,   we  obtain 


i^ 


^^ 


n+i 


n+i 


=     h 


H  h  ^'' 


1=0 


\p 


u 


n+i 


^n+i 


i^ 


=     hXI^i 


A^      0 


0         A, 


U   ^.     il 
n+i 


V. 


n+i 


=     hJZ^± 


~0 


A,    U   ^. 
1     n+i 


^2  ^n+i    i 


Recalling  that 
A, 


a 


11 


h^ 


(as  h  -*■  0) 


and 


,      det  A   /    i_   ^  N 

A^  'v  (as  h  -►  0) 

2      a^^ 


and  a, -J  <  0,   we  can  conclude  that  M--stability  is  a  suffi- 
cient condition  to  ensure   lu  I  -^  0  and   Iv  \      remains 

'  n'  '  n' 
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bounded   as      h  -^  0,      x      =  nh .      Since      P      is   a  "bounded  matrix 

n 

for  h  sufficiently  small,   |Y  |   and   |Z  |   will  remain 
■bo\inded. 

THEOREM  6.1.   Let  the  consistent  finite  difference 
equation 


(6.13)  ZZ  «i 


1^ 


^n+i 


^n+i 


i=0   ^   " 


n+i 


n+i 


«n    ! 

+ 

4 

where 


A 


n 


^11    ^12 


h^ 


h^ 


^21    ^22 


a  constaxit  matrix,  satisfy  the  following  conditions: 

a)  R^  =   (^(h^^^"^^)   ajid  R^   =   ^(h^"^)   as  h--0. 

b)  The  finite  difference  equation  is  M--sta.ble. 


c)  a.^^     <      0 


d)   Y^  -  y(x^) 


as  h  -+  0 


1-7 
i  =  0,1, .. . ,k-l 


{^(h^"'^)   and  Z.  -  z(x.  )   -  ^  (h"^"^) 


Then  the  difference  equation  is  convergent. 

Proof:   The  difference  equation  satisfying  the  error 


terms  e^  =  Y^  -  y(x^)   and  d^i  "  ^n  ~  ^^^n^   ^^ 
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(6.14)   XI  a. 


1^ 

3^0 


n+1 


n+i 


k 


where  the  last  term  is 
Introducing 


h  /^  ^i  \ 
i=0  ^  ^ 


'n+i 


d  ^.  " 
n+i  II 


+ 


;,  n  n  ,; 

il 

il  r2_^2  il 

ii  n  n  I, 


as   h  — ►  0. 


n 


n 


=  P 


u 


n 


V 


n 


where  P  was  defined  earlier  and 


P-^A^P   = 


A 


1    0 
0    A. 


where 


■^  h 


as  h  -^  0 


and 


■K        ^   iet_A,   ^g   ^  ^  0 


a 


11 


we  can  write  (6.l4)  as 


(6.15)  (cx^l  -  hP^ 


\  ° 

1 
) 

%+i 

0      ^2 

^n+i 

=  p 


■1 


R- 


n 

2 

R   -  T, 


n    n 


123 


where  we  have  used  the  explicit  form  of   P" 


0(1)  ^1) 
0{h'^)   ^(1) 


The  equation  for  u   can  "be  treated  like  the  one 
in  Theorem  2.2,  while  the  equation  for  v   follows  the 
standard  ODE  theory.   Q.E.D. 

THEOREM  6.2.   Let  the  consistent  finite  difference 
equation 


k 


(6.16) 


i=0 


a. 
1 


n+i 


n+i 


=  h 


P. 


i-0 


h"^F 


n+i 


n+i 


+ 


where   F  =  f(Y  ,Z  ,x  )   and   G  =  g(Y  ,Z  ,x  )   and  R 
n    '  n'  n'  n        n   ^^  n  n  n        n 

2 
and  R   are  the  round-off  errors  satisfy  the  following 

conditions: 

a)  R^     =     ^  (h^"^^)   and  R^   -  (f- {h^~^)      as   h  ^  0. 

b)  The  finite  difference  equation  is  relatively  stable. 

c)  The  first  and  second  partial  derivatives  of  f  and  g 
are  continuous  and  bounded  for  0  ^  x  ;^  1,   -oo  <  y  <  oo, 

and  -00  <  z  <  00.   Furthermore,   -L  <^  --—  <_  -L  <  0  in  this 
strip. 

d)  e.   =   {^(h-^"^)   and  d.   =  0  {h^~'^)      as  h  -►  0, 


\)  y    J~  y     •  •  •    f   ri'~  J-  • 
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Then  there  exists  a  constant   C   for  any  h   such  that 

|e  I  <  C  h^"^  and   |d  |  <  C  h^'^  for  n  =  0,1,- ••,N 
'  n'  —  '  n'  —  >    )         ) 

where  x„  =  1. 


Proof:   The  exact  solution  of  the  system  (6.2)  satisfies 


(  6 .  17 )  Z_  ^ 


y, 


n+i 


n+i 


n+i 
^+i 


+ 


where  f^+i  =  ^^^n+i' ^n+i'^n+i^  ^^^  S^+i  =  S^^n+i' ^n+i' 
^n+i^'  T^  =  (:>(h^"^^)  and  T^  =  0(h^~^)  as  h  --  0  by 
consistency  and  condition  c). 

Subtracting  (6.17)  from  (6.l6)  and  letting  ®n  "^  ^n  "  -^n 
and  <in  ^  ^n  '  ^n  Sives 


k 

>   "i 

^n-fi 

i=0     -^ 

d      . 

"n+i 

=  h 


1=^. 


i=0 


^  n+i  n+i 


n+i  °n+i 


r1  -  T"^ 

n  n 

2  2 

n  n 


(6.18) 


=  h 


P. 


i-0 


ry^ 


Yi   'f      ,  .  h  ^f 


7-. 


■yn+i 


zn+i 


g 


yn+i 


g 


zn+i 


n+i 


n+i 


+ 


r1  _  T-^ 

n  n 

2  2 

R^  -  T"^ 

n  n 


l-y 
T    ...  ,T    -2n  h  ^ 

Letting  N^  =  in  r 


we  can  obtain  a  bound  for  the 


error  terms  in  the  standard  way  (see  Theorems  4.1  and  5-2) 


for  0  <  n  <  N  .   Assuming  this  has  been  done  now  add 
—   —  o 


-.^-^  y:  Pi 


i=0 


af(y(x^),z(x^),x^) 


o 


o' ''  o'    e  ,  .  =  h 
n+i 


i^ 


'-'  H  Pi 


L  e  ^. 
o  n+i 
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to  both  sides  of  (6.18). 

Note:   If  •^—  is  a  monotonia    increasing  function 
ay  _  ^ 

as  is  usually  the  case,  then  L  =  L  . 


(6.19) 


i~0 


y~  (a.    +  h-'-'^p.L    )e   ^. 

4—yr    ^      X  1     O        n+1 


k 


i=0 


ex.    d    ,  . 
1     n+i 


h 


i=0 


h"^(L  +f      , . )    h  ^f      ^. 
^    o     yn+i  zn+i 


'yn+i 


g 


zn+i 


n+i 


^n+i 


+ 


1  ml 


R"    -    T 


R 


n  n 


1 
% 

2 


The  solutions  to  (6.19)  are 


r  n-k 
4 — TT  -^n-i-l  ^1    n 


(6.20) 


e   =  { 

n 


i=0 

,1 


VJ^- 


n 


n  >  k 


n  <  k 


where  $   and  t       are  defined  like  ^   and  ^   in 
Theorem  k .1  and 
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(6.21) 


(n-K 


1=0 


n 


T 


n 


n  <  k 


where  ¥   and  ^   are  defined  like  g   and  0   In 
■'ki                   n  °n       n 

Theorem  2.2.   By  continuity  and  condition  c)  for  some 
e  >  0 


o    yn+i '   —      ^   n+i    n+i 


for  0  <  n  <_  IT,  -  k  where  E  ^   max   |  e .  |    and 

^  ^       1=0,1, .^.,n 

D  =  max   I  d .  I  . 
^       1=0,1, .^.,n 

Using  (6.20)  and  following  the  procedure  of  Theorem 
^.1  we  can  obtain  a  bound  for   |s  I • 


ej   1   BK^E^  +  K^D^  +  (t  (  (E^  4-  D^Ej  +  ^  (h^'^') 


for  N  <  n  <  N^ 
o  —   —  1 


or  If  eK  <  1,   then 


(6.22)   E^  1  K^D^  +  (>((E^  +  D^)E^)  +  (^{h^'^) 


for   0  1  n  <_  N, 


where  K,,  K^,      and  K,   are  positive  constants 
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Letting  max  IJ^ 1  =     '^ >      we  have  from  (6.21)  that 
n-k   k 


n-k 

H 

1=0 


+  I^  ZZ  ^(h^"^)  +  l^nl 


k 

-r-P  "^ 
<    X 


t  TZ    iPjIt'^l^n  ■"  ^2°n^  +^(h^-^) 


for   0  ^  n  <  N-^^ 


where   |g   1  <  G^  ,   Ig^^l   1  ^2'   ^^^  x^  =  nh  <  x. 
Substituting   (6.22)  into  (6.23)  gives 

Choosing  x  small  enough  so  that  1  >  x  K^  we  can  write 

(6.24)  D^  1  (^(h^"^)  +  ^((E^  +  D^)En)    0  1  n  1  N^ 
and  consequently  from  (6.22) 

(6.25)  E^  1  (V(h^'^)  +  ^((E^  +  D^)E^)  . 

Adding  (6.2^1)  and  (6.25)  and  letting  W^  =  E^  +  D^ 
we  obtain 
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\     1     ^(h^"^)   +  ^(W^) 


The  proof  is  completed  like  Theorem  ^.1 
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VII .   APPLICATIONS  AND  PINAL  STATEMENTS 


A  series  of  non-linear  boundary  layer  problems  of  the 
kind  discussed  in  Chapter  4  were  solved  on  the  IBM  709^ 
and  the  CDC  6600  computers  by  various  finite  difference 
schemes.   At  first   7  was  chosen  equal  to  •^  and  mesh- 
sizes  of  h  =.01,  .001,  .0005,  .0005,  .0001,  and  .00005- 
Although  various  different  schemes  were  used  to  solve  the 
same  boundary  layer  equation,  we  concentrated  on  finite 
difference  methods  that  illustrated  the  need  for  the  new 
kind  of  stability  called  ix-stability,  i.e.,  most  of  our 
schemes  had  two  or  more  roots  of  p(C)   located  on  the 
unit  circle. 

We  also  considered  the  cases  of  double  roots  of  p(C) 
and  other  methods  that  were  stable  but  not  ix-stable.   The 
latter  methods  produced  divergent  results,  but  the  diver- 
gence was  not  clearly  obvious  \intil  x   was  outside  the 
boundary  layer.   Inside  the  boundary  layer  we  observed  two 
or  three  place  accuracy  with  the  exact  solution  or  with  the 
comparison  with  optimal  methods  that  were  ii-stable  schemes. 

Finally,  we  allowed   7  to  vary  between  0  and  1  and 
observed  the  effect  on  convergent  of  M.-stable  consistent 
schemes.   We  tried  7  =  .1,  .25,  .5,  .75.  and  .9-   We  noted 
that  for   .5  <  7  <  1  the  iterative  methods  required  much 
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more  time  ajid  a  few  cases  even  required  choosing  smaller 
values  of  h.   However,  consistent  ix-stable  methods  were 
observed  to  converge  even  for  these  values  of  y. 

The  accuracy  of  the  results  varied  with  the  precision 
of  the  difference  methods.  We  fo\and,  for  example,  that 

Euler's  method  (Y^^.]_  "  ^^  =  ^'"^   ^^^n'^n^   ^^"^®  °^-^"^ 
three  to  four  place  accuracy  while  optimal  methods  gave 
six  or  greater  place  accuracy  when  compared  with  known 
solutions . 

In  every  case  when  h  became  sufficiently  small 
schemes  which  were  predicted  to  converge  by  the  theory 
did  so,  while  schemes  which  were  predicted  to  diverge  over- 
flowed in  the  computer. 

We  illustrate  some  of  the  numerical  results  in  the 
following  graphs  smd  tables.   Graph  1  shows  the  numerical 
calculations  in  the  boundary  layer  for  a  particular  non- 
linear differential  equation  with  three  different  choices 
for  |x. 

Graph  2  compares  the  numerical  data  from  a  pL-stable 
method  with  a  stable  method.  In  this  case  the  M'-stablility 
condition  is  that  p,  <  1.  Notice  that  divergence  doesn't 
occur  with  the  stable  (but  not  M-'Stable)  method  inside  the 
boundary  layer,  but  outside  of  it. 

Figures  1,  2,  and  5  give  the  n\amber  of  exact  signifi- 
cant figures  in  the  numerical  data  when  compared  to  the 
exact  answer  calculated  to  eight  significant  figures.  The 
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methods  refer  to  the  following  numerical  schemes: 


Method  1   (Euler's  Method) 


^n+l  -  ^n  =  h'^^  P(Y^'^n) 


Method  2   (Adam's  Method) 

,1/2 
V2  -  Vl  =  V(5F(Y^^2'^n+2)  "^  ^^^  Vr^n+1^  "  ^^^n'^n^^ 


Method  3      (ii-stable  but  not  relatively  stable) 


Method  4      (all     p.'s   >  0;      roots      C   =  1,    -.2) 


Method  5      (A  poor     |JL-stable  method) 

\-i-2   -  ^n     =     "^"-^^   (  -5   ^^+2   +  ^n+1  "^   '^   ^n^ 


Method  6  (Simpson's  Rule;  divergent) 


This  data  was  calculated  on  the  IBM  709^,  and  it  was 
not  as  good  as  corresponding  data  calculated  on  the  CDC  6600. 
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We  suspect  that  round-off  errors  in  the  IBM  machine  caused 
the  poorer  showing. 

Figure  '4   compares  Method  1  (Euler's  Method)  and 
Method  7--the  optimal  method  for  the  roots   C  =  1^  -5,    -5: 


V3  -  2  V2  +  1-25  Y^+i  -  -25  Y^ 


h 


1/2 


(  ^1  F.^-.  +  19  P.^o  -  53  F„^,  +  17  F^) 


■^^  ^  "-^  ^n+3  ■"  ^^   ^n+2  -  ^^  ^n+1 


n 


This  data  was  calculated  on  the  CDC  6600. 

Finally,  in  Figure  5  we  compare  the  "worst  methods" 

\ — ^ 
(large  |C^|  or  large  2 i^o'^i'^  ^^^^   *^®  "best  methods" 

(those  guaranteed  to  converge  by  Corollary  1  or  Corollary  2 
of  Theorem  ^ .1) . 

k 
Method 


8  (large  XI  l^i  I  ) 


i=0 
V3  -  ^n+2  =  ^'^'  ^Fn+3  ^  ^°  ^n4-2  +  ^  F^^^  -  12  F^) 


Method  9   (large   [C^l;   roots   1,  -9,  -9  ) 
Vj  -  2-8  Y„^2  +  2.61  Y^^,  -  .81  Y^   = 

Methods  10  to  19 

The  optimal  methods  for  the  roots  C  =  1>  -•!,  -•!  and 
i  =  0,1,. ..,9. 
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For   i  =  0  and   1  Corollary  2  can  be  used  to  guarantee 
convergence . 

For   1=2,3,..., 9   all  p.'s   >   0.   Hence,  by  Corollary  1 
convergence  is  assured. 
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CHART  OP  SIGNIFICANT  FIGURES 
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CHART  OP  SIGNIFICANT  FIGURES 
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CHART  OF  SIGNIFICANT  FIGURES 
f(y,x)   =  -y(y  -  1)(20  x  +  10);   y(0)  =  2 


Method 

x-coordinate 

8 

9 

10-19 

Comment 

.0001 

1 

1 

9 

.0005 

2 

1 

9 

.001 

2 

1 

9 

.005 

4 

3 

9 

.01 

6 

3 

9 

.015 

2 

3 

9 

End  Boimdary  L&yi 

.02 

3 

3 

9 

.03 

0 

5 

9 

Divergence  in  8 

.1 

0 

1 

9 

.2 

0 

9 

Divergence  in  9 

.3 

0 

9 

A 

9 

•5 

9 

.6 

9 

•7 

9 

.8 

9 

.9 

9 

1.0 

9 

Figure  5 


Data  from  CDC  6600 


139 


BIBLIOGRAPHY 


[1]   Dahlquist  (1956):   Convergence  and  stability  in  the 

niomerical  integration  of  ordinary  differential  equations 
Math.  Scan.,  ^4,   33-55- 

[2]   Henrici  (1962):   Discrete  Variable  Methods  in  Ordinary 
Differential  Equations,  John  Wiley  and  Sons,  New  York. 

[3]   Hull  and  Luxemburg  (1960):   Numerical  methods  and  exis- 
tence theorems  for  ordinary  differential  equations. 
NiAm.  Math.,  2,  30-41. 

[k]      Vasil'eva  (1963):   Asymptotic  behavior  of  solutions. 
Russian  Math.  Surveys,  I8,  13-84. 

[5]   Hull  and  Newbery  (1959):   Error  bounds  for  a  family 
of  three-point  integration  procedures.   J-  Soc . 
Indust.  Appl.  Math.,  7,    402-4l2. 

[6]  Hull  and  Newbery  (1961):  Procedures  which  minimize 
propagated  errors.  J.  Soc.  Indust.  Appl.  Math.,  9, 
28-47. 


140 


Graph 


Behavior  of  Solutions  in  the  Boundary  Layer 


)uy'  =  -y(y-i)(20x4-  lO) 
y(0)=2      hA=h'/2 


o 

Q. 

£ 
o 
o 


141 


ro 


Cl 

o 


V) 

E 

JZ 

o 

(/) 

Qi 

Xi 
D 

I 

•o 

c 
o 

0) 

o 
en 


c 
o 

V) 


o 

o 

o 
o 

+ 

X 

II 

-C 

o 

OJ 

OJ 

-C 
II 

1 

i. 

1 

II 

o 

f 

\ 


s 


CO 

+ 

00 

(VJ 

SI 

II 
c 
> 
I 
CJ 

+ 


Si 

—  JD 


^  00 
Q  ^ 

+    + 


o 

c 

o 
CO 


^ 


of 
in 


4- 
00 

cr> 

+     _a; 
o 


in 


(\j 


II 

c 

>- 

I 
(\J 


CO 
I 


(\j 


o 


CD 


CD 


C\J 


CD 


CD 


l42 


This  report  was  prepared  as  an  account  of 
Government  sponsored  work.   Neither  the 
United  States,  nor  the  Commlaslon,  nor  any 
person  acting  on  behalf  of  the  Connnlsslon: 

A.  Makes  any  warranty  or  representation, 
express  or  Implied,  with  respect  to  the 
accuracy,  completeness,  or  usefulness  of 
the  Information  contained  In  this  report, 
or  that  the  use  of  any  Information, 
apparatus,  method,  or  process  disclosed 
in  this  report  may  not  Infringe  privately 
owned  rights;  or 

B.  Assumes  fiuiy  liabilities  with  respect  to 
the  use  of,  or  for  damages  resulting  from 
the  use  of  any  Information,  apparatus, 
method,  or  process  disclosed  In  this 
report. 

As  used  In  the  above,  "person  acting  on  behalf 
of  the  Commission"  Includes  any  employee  or 
contractor  of  the  Commission,  or  employee  of 
such  contractor,  to  the  extent  that  such  em- 
ployee or  contractor  of  the  Commission,  or 
employee  of  auch  contractor  prepares,  dis- 
seminates, or  provides  access  to,  any  Infor- 
mation pursuant  to  his  employment  or  contract 
with  the  Commission,  or  his  employment  with 
such  contractor. 
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